The myometrial transcriptome changes in mares with endometrosis

Sequencing of cDNA libraries constructed from total RNA isolated from the myometrium of mares with different endometrial scores (categories I, IIA, IIB, and III) provided 383.8 million raw reads ranging from 18.1 to 20.6 million per sample. After removing low-quality reads (reads length < 50 bp, Phred score Q < 30), the remaining 364.2 million high-quality reads (from 17.2 to 19.6 million per sample) were mapped to the horse reference genome. The number of reads aligned uniquely to the reference genome ranged from 16.2 to 18.5 million per sample, and an average of 91.1% of these reads were mapped to unique locations. The total number of expressed genes in the mares’ myometrium of all examined samples ranged from 15,995 to 20,208 (Supplementary Table 1). The results of the PCA and sample distance revealed that samples were formed as separate groups (Fig.  1 A,B). The Volcano plot presents the significant differences (P-adjusted < 0.05, |log2FC|≥ 1.0) in expression profiles of transcripts that altered in the myometrium of mares with different stages of endometrosis (IIA vs I, Fig.  1 C; IIB vs I, Fig.  1 D; III vs I, Fig.  1 E; IIB vs IIA, Supplementary Fig. 1A; III vs IIB, Supplementary Fig. 1B; III vs IIA, Supplementary Fig. 1C).The most marked differences in gene expression patterns were found between the myometrium of mares with endometrium category I compared to categories IIA, IIB, and III; Fig.  1 ). Fig. 1 Global changes in transcriptional profiles of mares with different stages of endometrosis. ( A ) Graphical presentation of the first (PC1) and second (PC2) principal components of the samples. ( B ) Samples distance matrix for 500 most variable genes. ( C ) The volcano plot presenting differentially expressed transcript active regions (deTARs; P-adjusted < 0.05. and |log2FC|≥ 1) in mares with different stages of endometrosis IIA vs I, ( D ) IIB vs I, and E ) III vs I. TARs are represented by multicolored circles, where red color represents up − regulated deTARs and green down − regulated deTARs. Gray circles represent all genes identified in the myometrium with no significant changes and red/green triangles deTARs with -log10(P-adjusted) > 80 to increase readability. Global changes in transcriptional profiles of mares with different stages of endometrosis. ( A ) Graphical presentation of the first (PC1) and second (PC2) principal components of the samples. ( B ) Samples distance matrix for 500 most variable genes. ( C ) The volcano plot presenting differentially expressed transcript active regions (deTARs; P-adjusted < 0.05. and |log2FC|≥ 1) in mares with different stages of endometrosis IIA vs I, ( D ) IIB vs I, and E ) III vs I. TARs are represented by multicolored circles, where red color represents up − regulated deTARs and green down − regulated deTARs. Gray circles represent all genes identified in the myometrium with no significant changes and red/green triangles deTARs with -log10(P-adjusted) > 80 to increase readability. In the myometrium of mares with mild endometrosis (IIA vs I), 665 deTARs were evaluated [469 up-regulated, including 361 protein-coding deTARs, i.e., differentially expressed genes, DEGs, and 196 down-regulated, including 180 DEGs; Table 1 , Supplementary Table 2 (IIA vs I)]. In mares with moderate endometrosis (IIB vs I), 491 deTARs [328 up-regulated, including 251 DEGs, and 163 down-regulated, including 119 DEGs; Table 1 , Supplementary Table 3 (IIB vs I)] were identified. In mares with severe endometrosis (III vs I), 499 deTARs [338 up-regulated, including 267 DEGs, and 161 down-regulated, including 144 DEGs; Table 1 , Supplementary Table 4 (III vs I)] were found. In the myometrium of mares in each stage of endometrosis (IIA vs I, IIB vs I, and III vs I), 200 common deTARs were identified [DEGs: 131 up-regulated, 21 down-regulated; lncRNAs: 46 up-regulated, and one down-regulated; one up-regulated pseudogene; Table 1 , Supplementary Table 5 (Common deTARs)]. Figure  2 shows a graphical representation of the number of evaluated deTARs in each comparison. The expression profile of the top 50 deTARs (i.e., deTARs with the lowest P-adjusted value) for IIA vs I, IIB vs I, III vs I, and all common deTARs between comparisons is presented in Fig.  3 . The top 20 down-regulated and top 20 up-regulated (P-adjusted < 0.05, and |log2FC| ≥ 1) common deTARs identified in the myometrium of mares with endometrium categories IIA, IIB, and III vs I are presented in Table 2 . Table 1 Differentially expressed transcriptionally active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) identified in the myometrium of mares during the mid-luteal phase of the estrous cycle with different endometrial score—summary results. Condition Biotype Total Up-regulated Down-regulated IIA vs I IG_V_gene 1 1 0 TR_C_gene 1 1 0 lncRNA 119 104 15 miRNA 1 1 0 protein_coding 541 361 180 pseudogene 1 1 0 snoRNA 1 0 1 SUMMARY 665 469 196 IIB vs I lncRNA 116 74 42 miRNA 1 1 0 protein_coding 370 251 119 pseudogene 2 1 1 snRNA 1 1 0 snoRNA 1 0 1 SUMMARY 491 328 163 IIB vs IIA lncRNA 37 6 31 protein_coding 92 46 46 snoRNA 1 0 1 SUMMARY 130 52 78 III vs I lncRNA 86 69 17 protein_coding 411 267 144 pseudogene 1 1 0 snRNA 1 1 0 SUMMARY 499 338 161 III vs IIA lncRNA 11 6 5 protein_coding 31 19 12 snoRNA 1 1 0 SUMMARY 43 26 17 III vs IIB lncRNA 44 39 5 miRNA 2 0 2 protein_coding 62 37 25 snoRNA 2 2 0 SUMMARY 110 78 32 Common IIA, IIB, and III vs I lncRNA 47 46 1 protein_coding 152 131 21 pseudogene 1 1 0 SUMMARY 200 178 22 IIA vs I (without common) IG_V_gene 1 1 0 TR_C_gene 1 1 0 lncRNA 72 58 14 miRNA 1 1 0 protein_coding 389 230 159 snoRNA 1 0 1 SUMMARY 465 291 174 IIB vs I (without common) lncRNA 69 28 41 miRNA 1 1 0 protein_coding 218 120 98 pseudogene 1 0 1 snRNA 1 1 0 snoRNA 1 0 1 SUMMARY 291 150 141 III vs I (without common) lncRNA 39 23 16 protein_coding 259 136 123 snRNA 1 1 0 SUMMARY 299 160 139 Fig. 2 The differentially expressed transcript active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) identified in the myometrium of mares with different stages of endometriosis: ( A ) IIA vs I, IIB vs I, and III vs I pairs; ( B ) IIB vs IIA, III vs IIB, and III vs IIA pairs. Fig. 3 Heatmap of the top 50 identified differentially expressed transcript active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) in the myometrium of mares with different stages of endometrosis: ( A ) IIA vs I, ( B ) IIB vs I, ( C ) III vs I, ( D ) Common deTARs between comparisons (IIA vs I, IIB vs I and III vs I). The top genes from panels A-C were selected based on the lowest P-adjusted value. Table 2 The top 20 down-regulated and top 20 up-regulated among commonly differentially expressed transcriptionally active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) identified in the myometrium of mares during the mid-luteal phase of the estrous cycle with endometrium category IIA, IIB and III vs I. Ensembl Symbol Name Biotype Fold change IIA vs I IIB vs I III vs I Down-regulated TARs ENSECAG00000007493 ACAN aggrecan protein_coding -19.90 -18.14 -13.19 ENSECAG00000004696 ENSECAG00000004696 neuroendocrine secretory protein 55 protein_coding -17.47 -17.33 -12.05 ENSECAG00000058396 ENSECAG00000058396 NA protein_coding -8.27 -9.39 -9.16 ENSECAG00000046538 ENSECAG00000046538 NA lncRNA -5.22 -3.94 -3.28 ENSECAG00000013863 LEPR leptin receptor protein_coding -4.27 -2.34 -2.33 ENSECAG00000028239 ENSECAG00000028239 mammaglobin-A protein_coding -3.50 -3.86 -3.52 ENSECAG00000011663 NPTX1 neuronal pentraxin 1 protein_coding -3.33 -2.96 -4.68 ENSECAG00000010475 IL33 interleukin 33 protein_coding -1.83 -3.64 -2.15 ENSECAG00000008004 PDE10A phosphodiesterase 10A protein_coding -1.71 -2.19 -1.72 ENSECAG00000035602 ENSECAG00000035602 NA protein_coding -1.68 -2.98 -4.05 ENSECAG00000009311 SH3TC2 SH3 domain and tetratricopeptide repeats 2 protein_coding -1.65 -1.62 -1.35 ENSECAG00000000046 FAM13C family with sequence similarity 13 member C protein_coding -1.62 -1.95 -2.81 ENSECAG00000014487 SLC14A1 solute carrier family 14 member 1 protein_coding -1.57 -1.22 -2.91 ENSECAG00000023668 RGS5 regulator of G protein signaling 5 protein_coding -1.57 -1.43 -1.83 ENSECAG00000017227 GASK1B golgi associated kinase 1B protein_coding -1.53 -1.92 -1.81 ENSECAG00000012752 EPHA4 EPH receptor A4 protein_coding -1.50 -1.01 -1.64 ENSECAG00000008923 THBS1 thrombospondin 1 protein_coding -1.45 -1.40 -1.73 ENSECAG00000003655 ENSECAG00000003655 NA protein_coding -1.24 -2.29 -2.46 ENSECAG00000036379 PLN phospholamban protein_coding -1.17 -1.88 -2.94 ENSECAG00000018511 SORBS1 sorbin and SH3 domain containing 1 protein_coding -1.12 -1.49 -1.63 Up-regulated TARs ENSECAG00000048284 ENSECAG00000048284 NA lncRNA 22.41 20.83 25.65 ENSECAG00000049934 ENSECAG00000049934 NA lncRNA 22.42 22.12 19.40 ENSECAG00000041200 ENSECAG00000041200 NA lncRNA 22.43 29.20 32.09 ENSECAG00000055055 ENSECAG00000055055 NA lncRNA 22.51 24.13 20.95 ENSECAG00000022644 ENSECAG00000022644 C-type lectin-like domain family 1 protein_coding 22.52 26.17 22.32 ENSECAG00000051478 ENSECAG00000051478 NA lncRNA 22.65 15.99 23.68 ENSECAG00000018390 FFAR4 free fatty acid receptor 4 protein_coding 22.75 25.30 26.93 ENSECAG00000030129 C2orf16 chromosome 2 open reading frame 16 protein_coding 23.84 25.08 20.13 ENSECAG00000018546 TP63 tumor protein p63 protein_coding 24.68 24.50 19.27 ENSECAG00000037104 ENSECAG00000037104 NA lncRNA 25.10 24.93 24.93 ENSECAG00000053211 ENSECAG00000053211 NA lncRNA 25.53 26.97 26.65 ENSECAG00000045721 ENSECAG00000045721 NA lncRNA 26.50 24.97 28.53 ENSECAG00000058416 ENSECAG00000058416 NA lncRNA 26.78 32.23 27.75 ENSECAG00000035475 ENSECAG00000035475 N-acetyltransferase 8B protein_coding 28.32 4.64 28.05 ENSECAG00000028354 ENSECAG00000028354 NA lncRNA 30.10 17.08 37.10 ENSECAG00000024782 HYAL4 hyaluronidase 4 protein_coding 30.29 29.05 35.53 ENSECAG00000047247 ENSECAG00000047247 NA lncRNA 32.75 33.96 31.99 ENSECAG00000015436 PLEKHS1 pleckstrin homology domain containing S1 protein_coding 34.16 32.27 36.38 ENSECAG00000022699 PKHD1 PKHD1 ciliary IPT domain containing fibrocystin/polyductin protein_coding 37.90 37.71 39.06 ENSECAG00000047676 ENSECAG00000047676 NA lncRNA 46.69 42.56 41.92 Differentially expressed transcriptionally active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) identified in the myometrium of mares during the mid-luteal phase of the estrous cycle with different endometrial score—summary results. The differentially expressed transcript active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) identified in the myometrium of mares with different stages of endometriosis: ( A ) IIA vs I, IIB vs I, and III vs I pairs; ( B ) IIB vs IIA, III vs IIB, and III vs IIA pairs. Heatmap of the top 50 identified differentially expressed transcript active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) in the myometrium of mares with different stages of endometrosis: ( A ) IIA vs I, ( B ) IIB vs I, ( C ) III vs I, ( D ) Common deTARs between comparisons (IIA vs I, IIB vs I and III vs I). The top genes from panels A-C were selected based on the lowest P-adjusted value. The top 20 down-regulated and top 20 up-regulated among commonly differentially expressed transcriptionally active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) identified in the myometrium of mares during the mid-luteal phase of the estrous cycle with endometrium category IIA, IIB and III vs I. In the myometrium of mares with mild endometrosis (endometrium IIA vs I), 313 (190 up-regulated and 123 down-regulated) changes in gene expression were identified that do not appear in comparisons IIB vs I and III vs I [Supplementary Table 6 (exclusive IIA vs I), Table 1 , Fig.  2 A]. In the myometrium of mares with moderate endometrosis (endometrium category IIB vs I), there were identified 169 (69 up-regulated and 100 down-regulated) changes in gene expression that do not appear in comparisons IIA vs I and III vs I [Supplementary Table 7 (exclusive IIB vs I), Table 1 , Fig.  2 A]. In the myometrium of mares with severe endometrosis (endometrium category III vs. I), 161 (77 up-regulated and 84 down-regulated) changes in gene expression that do not appear in comparisons IIA vs I and IIB vs I [Supplementary Table 8 (exclusive deTARs III vs I, Table 1 , Fig.  2 A)]. In the myometrium of mares, 130 deTARs were identified in those with endometrium IIB compared to IIA (52 up-regulated and 78 down-regulated; Supplementary Table 9); IIB vs IIA, (Table 1 ), 43 deTARs were identified in those with endometrium III compared to IIA [(26 up-regulated and 17 down-regulated; see Supplementary Table 10 (III vs IIA), Table 1 )], and 110 deTARs were identified in those with endometrium III compared to IIB [(78 up-regulated and 32 down-regulated; Supplementary Table 11 (III vs IIA), Table 1 )]. Figure  2 B presents summary results of the number of deTARs identified between comparisons IIB vs IIA, III vs IIA, and III vs IIB. All evaluated GO terms for identified deTARs (up- and down-regulated together and separately in the myometrium with mild, moderate, and severe endometrosis) are presented in Supplementary Table 12 (GO terms). Supplementary Fig. 3 shows a graphical presentation of the number of evaluated GO terms for each comparison with an indication of the most enriched GO terms within BP, CC, and MF for up-regulated and down-regulated deTARs in the myometrium of mares with endometrium category IIA vs I, IIB vs I, and III vs I. The most important contribution of common DEGs for top identified GO terms is presented in Fig.  4 A, and deTARs assigned to GOs connected with muscle structure and contractions are presented in Fig.  4 B. Fig. 4 Gene Ontology (GO) analysis of evaluated differentially expressed transcript active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) in the myometrium of mares with different stages of endometrosis. ( A ) Sankey plot illustrating the most important contribution of common deTARs for top identified GO terms. ( B ) The selected GO terms involved in muscle functionality identified for up- and down-regulated deTARs in IIA vs I comparison, down-regulated deTARs in IIA vs I comparison, and down-regulated deTARs in III vs I comparison. Gene Ontology (GO) analysis of evaluated differentially expressed transcript active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) in the myometrium of mares with different stages of endometrosis. ( A ) Sankey plot illustrating the most important contribution of common deTARs for top identified GO terms. ( B ) The selected GO terms involved in muscle functionality identified for up- and down-regulated deTARs in IIA vs I comparison, down-regulated deTARs in IIA vs I comparison, and down-regulated deTARs in III vs I comparison. In brief, in the myometrium of mares with endometrium category IIA vs I up-regulated TARs were grouped into 75 GO biological processes (BP) terms, 11 GO cellular components (CC) terms, and 14 GO molecular function (MF) terms, and down-regulated TARs were grouped to 205 GO BP terms. In the myometrium of mares with endometrium category IIB vs I up-regulated TARs were grouped into 23 GO BP terms, six CC terms, and 19 GO MF terms, and down-regulated TARs were grouped into eight GO BP terms, one GO CC term, and two GO MF terms. In the myometrium of mares with endometrium category III vs I, up-regulated TARs were grouped into 41 GO BP terms, 14 GO CC terms, and nine GO MF terms, and down-regulated TARs were grouped into 160 GO BP terms and five GO CC terms. Among common deTARs (up- and down-regulated together) in the myometrium of mares with endometrosis (endometrium category IIA, IIB, and III) compared to mares without fibrotic changes in the endometrium (category I) were evaluated 32 GO BP terms, 15 GO CC terms and seven GO MF terms (Supplementary Table 13). The most enriched with evaluated deTARs among GO BP was adaptive immune responses (GO:0002250; 16 deTARs). Among GO CC the most enriched were cell surface (GO:0009986; 22 deTARs), supramolecular fiber (GO:0099512; 18 deTARs), and polymeric cytoskeletal fiber (GO:0099513; 17 deTARs). Among GO MF the most enriched one was oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen reduced flavin flavoprotein as one donor , and incorporation of one atom of oxygen (GO:0,016,712, 5 deTARs). All evaluated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways are presented in Supplementary Table 4 (all deTARs). TARs with lowered expression in the myometrium of mares with endometrium category IIA vs I and IIB vs I were categorized into eight and 12 KEGG pathways, respectively (Supplementary Table 14). In both comparisons, down-regulated TARs were categorized to cAMP signaling pathway (ecb04024; 10 deTARs in IIA vs I and eight deTARs in IIB vs I), focal adhesion (ecb04510; 10 deTARs in IIA vs I and seven deTARs in IIB vs I), calcium signaling pathway (ecb04020 10 deTARs in IIA vs I and seven deTARs in IIB vs I) and ECM-receptor interaction (ecb04512; seven deTARs in IIA vs I and five deTARs in IIB vs I) pathways (Fig.  5 ). Additionally, down-regulated TARs in IIB vs I were categorized to oxytocin signaling pathway (ecb04921; five deTARs) (Supplementary Table 14, Fig.  5 ). deTARs in category III vs I were categorized into serotonergic synapse KEGG pathway (ecb04726; four up-regulated TARs) and neuroactive ligand-receptor interaction KEGG pathway (ecb04080; nine down-regulated TARs) (Supplementary Table 14). Fig. 5 Selected KEGG pathways enriched by differentially expressed transcript active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) evaluated in the myometrium of mares with different stages of endometrosis. Selected KEGG pathways enriched by differentially expressed transcript active regions (deTARs; P-adjusted < 0.05, and |log2FC|≥ 1) evaluated in the myometrium of mares with different stages of endometrosis. The validation procedure confirmed the results obtained with RNA-seq. Changes in the expression of selected for validation genes determined with transcriptome profiling correspond with results obtained with qPCR. Figure  6 shows the relative expression of the selected validated genes and summarizes fold change values obtained after transcriptome profiling and qPCR for each validated gene. Fig. 6 Validation of RNA-seq results. Different lower-case letters (a b, c) indicate statistically significant differences (P < 0.05) in gene relative expression in the myometrium of mares with different endometrial scores. Data were analyzed using multi-way ANOVA with Fisher’s post-hoc test. ( A ) Summary results of identified changes ( B ) interleukin 33 ( IL33 ), ( C ) cellular communication network factor 1 ( CCN1 ), D) epidermal growth factor ( EGF ), E) collagen type IV alpha 1 ( COL4A1 ), and F) nectin-4 ( NECTIN4 ). Validation of RNA-seq results. Different lower-case letters (a b, c) indicate statistically significant differences (P < 0.05) in gene relative expression in the myometrium of mares with different endometrial scores. Data were analyzed using multi-way ANOVA with Fisher’s post-hoc test. ( A ) Summary results of identified changes ( B ) interleukin 33 ( IL33 ), ( C ) cellular communication network factor 1 ( CCN1 ), D) epidermal growth factor ( EGF ), E) collagen type IV alpha 1 ( COL4A1 ), and F) nectin-4 ( NECTIN4 ).

The Local Ethics Committee for Experiments on Animals in Olsztyn, Poland (Agreements No. 51/2011) approved all material collection procedures. The mares were euthanized in accordance with European Legislation (EFSA, AHAW/04–027) to eliminate pain and suffering. The experiments were performed on tissue samples collected post-mortem during the regular economical slaughter. All methods were carried out following relevant guidelines and regulations indicated where appropriate. The material was collected from clinically healthy Polish cold-blood mares (2–20 years) weighing 500 ± 100 kg, designated for regular economic slaughter at a local abattoir (Rawicz, Poland). The health of the animals was monitored by an official government veterinary inspector and by referral to historical health records for individual animals. The reproductive history of mares was not registered. The entire uteri with ovaries were collected post-mortem within 15 min of each mare’s death. For the study were selected mares with well-developed corpus luteum (CL) and accompanied by 15- to 20-mm diameter follicles, i.e. during the mid-luteal phase of the estrous cycle. Additionally, peripheral blood samples were collected into heparinized tubes immediately before the slaughter for further estrous cycle confirmation based on determination of serum progesterone (P 4 ) level by radioimmunoassay (RIA) (KIP1458, DIAsource, Belgium). Only samples that corresponded to peripheral blood serum P 4 level higher than 6 ng/mL were included in the study as originating from mares during the mid-luteal phase of the estrous cycle 56 , 57 . This phase of the estrous cycle was selected as the hormonal environment mimics the pre-implantation period when embryo migration occurs and disturbed contractile activity of the myometrium, caused by endometrosis, may contribute to early embryo death. The classification of mares according to the severity of inflammation and endometrosis, i.e. endometrium category I (without histopathological changes), IIA (mild condition), IIB (moderate condition), and III (severe condition) of the Kenney and Doig classification system 2 concerning inflammatory infiltrates, fibrosis, inflammatory infiltrates, and degree of dilatation of endometrial glands and lymphatic vessels was performed as reported previously 7 . Moreover, each specimen was evaluated for inflammation based on counting both mononuclear and polymorphonuclear neutrophil (PMN) cell infiltration of the endometrial luminal epithelium and stratum compactum, based on the criteria described by Ricketts, Ricketts and Alonso, and Nielsen 58 – 60 . Mares with endometritis were excluded from the study. Myometrial samples designated for RNA extraction were retrieved from uterine horns from the site of the active ovary. The perimetrium and the myometrium were excised, the cross-sections of the myometrium were picked and immediately placed in cryotubes containing RNAlater™ (#AM7021; Invitrogen, USA) and kept at a 4 °C overnight. RNAlater™ was removed from the samples 24 h after tissue immersion. The samples were then stored at -80 °C until use. All the analyses were performed on myometrial samples obtained from mares with endometrium category I (n = 5), category IIA (n = 5), category IIB (n = 4), and category III (n = 5). Total RNA was extracted from 20 mg fragments of myometrial tissue using TRI Reagent ® (T9424-200 ML, Sigma Aldrich, Germany) according to the manufacturer’s instructions after the tissue samples were homogenized on ice using a Tissue Ruptor homogenizer (Qiagen, USA). Aliquots were used for (1) initial RNA quality control (optical density, OD, A260/A280) and determination of RNA concentration using NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific, USA), (2) microfluidic electrophoresis for determination of RNA integrity number (RIN; 28 S/18 S ratio) using Bioanalyzer 2100 with RNA 6000 Nano LabChip kit (Agilent Technologies, USA), (3) construction of cDNA libraries and next-generation sequencing (RNA-seq) using NovaSeq6000 System (Illumina, USA), and (4) cDNA synthesis for further validation step using qPCR. Only samples that possessed OD 1.8–2.0 and RIN ≥ 7.5 were proceeded further. Construction of cDNA libraries and RNA-seq was done by an outsourcing company (Macrogen Europe, Netherlands). In brief, there were constructed TruSeq mRNA stranded libraries, and only libraries that achieved ≥ 3 nM were selected for sequencing. Run configurations for sequencing were 2 × 100 bp and a throughput of 40 M paired reads per sample. The raw data were submitted to the NCBI BioProject under accession No. PRJNA1023067, for further transcriptome profiling and bioinformatic analysis of gene expression. The quality of raw reads was evaluated using FASTQC ( https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ ). Low-quality reads and adapters were removed using Trimmomatic (version 0.39) using LEADING:6 TRAILING:6 SLIDINGWINDOW:5:20 MINLEN:50 AVGQUAL:30 parameters 61 . The obtained reads were mapped to the horse reference genome (EquCab3.0; Ensembl release 109) using STAR software (version 2.7.10a) 62 . After reads mapping, raw counts per gene were calculated by featureCounts (version 2.0.3). The differentially expressed transcript active regions (deTARs) as well as the corresponding P-adjusted values were determined employing R statistical software (version 4.2.1.) using the DESeq2 package (version 1.36.0) 63 . The evaluation of deTARs comprised of comparison of the myometrial transcriptomes of mares with mild, moderate, and severe endometrosis vs mares with healthy endometrium (IIA vs I, IIB vs I, and III vs I), and among endometrosis stages (IIB vs IIA, III vs IIA, III vs IIB). The Independent Hypothesis Weighting (IHW) package 64 was used to increase the statistical power while controlling the false discovery rate (FDR). Additionally, the log 2 fold change (log2FC) shrinkage was estimated by the ashr package 65 . The threshold for the significantly different expression was set at P-adjusted < 0.05 and log2FC ≥ 1.0 or log2FC ≤  − 1.0. The visual presentation of the results was performed by R software using ggplot2 (version 3.3.6; ISBN: 978–0-387–98,141-3) package. The heatmaps for gene expression were generated by the pheatmap package (version 1.0.12; Pheatmap: pretty heatmaps. R package version, 1(2), p.726). Clustering was done using Euclidean distance. The evaluated deTARs were divided according to their biotypes into protein-coding deTARs, i.e. differentially expressed genes (DEGs) and non-protein coding, uncovered regulatory region groups, including intraclonal diversification of immunoglobulin (Ig) variable (V) genes (IG_V_genes), Immunoglobulin (Ig) variable chain and T-cell receptor (TcR) genes (TR_C_genes), long non-coding RNAs (lncRNAs), microRNAs (miRNAs), small nuclear RNAs (snRNAs), small nucleolar RNAs (snoRNAs), and pseudogenes. Prior to functional analysis, the human gene orthologs were identified for equine DEGs utilizing biomaRt (version 2.52.0) 66 . Orthologs of genes were used for functional analysis based on the Gene Ontology (GO) database using clusterProfiler (version 4.4.4) 67 , DOSE (version 3.22.1) 68 , and org.Hs.eg.db (version 3.16.0) 69 packages of R software, with the established criterion P-adjusted < 0.05. Additionally, the Kyoto Encyclopedia of Genes and Genomes (KEGG database) was used to ascribe identified DEGs to particular biological mechanisms and cellular pathways (the established criteria: P-adjusted < 0.05) 70 – 72 . The KEGG enrichment analysis was performed by the clusterProfiler, DOSE, and org.Hs.eg.db packages of R software. The visual presentation of the results was performed by R software using ggplot2. To obtain cDNA, 1500 ng of total RNA aliquots were used. To eliminate potential gDNA contamination, each sample was treated with a DNA wipeout buffer containing Dnase I (DNase 1 Amplification Grade, AMPD1-1KT, #051M6157, Sigma Aldrich) following a manufacturer’s protocol. Afterward, cDNA was synthesized using a High Capacity cDNA Reverse Transcription Kit (#4,368,814, Applied Biosystems, USA), following a manufacturer’s protocol. Samples were stored at – 20 °C until further use. The amplification was conducted in a Viia7 system (Applied Biosystems). The reaction mixture for the qPCR assay consisted of 3 μL DNA (2.5 ng), 5 μL TaqMan Universal PCR Master Mix (#4,440,049, Applied Biosystem), 0.5 μL TaqMan probes (Thermo Fisher Scientific), and 1.5 μL nuclease-free water (#129,114, Qiagen, USA) to a final volume of 10 μL. For validation, protein-coding genes that were altered in the myometrium of mares in each stage of endometrium classification, i.e. common DEGs ( interleukin 33 , IL33 , assay ID Ec06977210_m1), or were altered in a specific stage of the condition, i.e. IIA vs I ( epidermal growth factor , EGF , assay ID Ec06967921_m1, and cellular communication network factor 1 , CCN1 , Ec07060486_g1), IIB vs I ( collagen type IV alpha 4 chain , COL4A4 , assay ID Ec06946014_m1) or III vs I ( nectin cell adhesion molecule 4, NECTIN4 , assay ID Ec07034730_m1) were selected. To select the best combination of reference genes the stability of glyceraldehyde-3-phosphate dehydrogenase ( GAPDH , assay ID: Ec03210916_gH), hypoxanthine–guanine phosphoribosyl transferase 1 ( HPRT1 , assay ID: Ec03470217), ubiquitin-conjugating enzyme E2 B ( UBE2B , assay ID: Ec04660756), ribosomal protein S18 ( RPS18 , assay ID: Ec06969343), succinate dehydrogenase complex, subunit A, flavoprotein variant ( SDHA , assay ID Ec03470487_m1), and ribosomal protein L32 ( RPL32 , assay ID: Ec06951800) were tested using the NormFinder software program (MOMA, Aarhus University Hospital, Denmark) as previously described by Andersen et al. 73 . SDHA , UBE2B , and RPL32 were chosen as the best combination of reference genes for myometrial tissue of mares during the mid-luteal phase of the estrous cycle. The stability value for RPL32 was 0.167, SDHA was 0.198, and UBE2B was 0.343. Amplification was held as follows: initial denaturation (10 min at 95 °C), 45 reaction cycles of denaturation for 15 s at 95 °C, and annealing for 1 min at 60 °C. Control reactions without a template were performed to confirm that products were free from genomic DNA contamination. The qPCR results were analyzed using the method described by Zhao and Fernald 74 . Gene expression data were expressed as relative to those of the best combination of reference genes and were presented as arbitrary units ± standard error of the mean (SEM).

The current study provides, for the first time, the transcriptomic characterization of the myometrium in mares with endometrosis for a better understanding of the molecular background of the histological and functional changes that occur in the tissue during this condition. We showed that the expression of 665 TARs was altered in the myometrium of mares with endometrium assigned to category IIA compared to category I, 491 deTARs in IIB compared to category I, and 499 deTARs in III compared to category I. Additionally, 200 deTARs were shared in the myometria of categories IIA, IIB, and III compared to category I. The majority of the observed common changes were identified as increased expression and corresponded to protein-coding deTARs (DEGs). The myometrial DEGs in category IIA vs I and IIB vs I were annotated to KEGG pathways including the calcium signaling pathway , cAMP signaling pathway , ECM-receptor interaction , and focal adhesion KEGG signaling pathways. Additionally, the myometrial DEGs in category IIB vs I were annotated in the oxytocin signaling pathway . Furthermore, it was demonstrated that especially in mares in the early stage of endometrosis myometrial DEGs with lowered expression were associated with numerous biological processes related to the contractile activity such as muscle contraction , regulation of muscle contraction , and regulation of muscle contraction GO BP terms, as well as muscle development and proliferation including muscle structure development , muscle tissue development , muscle cell proliferation , muscle cell differentiation , and actin filament organization GO BP terms. Myometrial contractile activity ensures the transportation of gametes, embryo migration, and accurate placement prior to implantation 10 – 14 , 19 . A proper myometrial function is also required for adequate uterine clearance after mating or labor, preventing an excessive inflammatory response in the uterus, which in mares can lead to the development of persistent endometritis and consequently, endometrosis 8 , 9 , 16 , 20 . As previously reported, histological changes in the myometrium and reduced contractile activity of the myometrial tissue have been observed in mares with endometrosis 15 , 21 . These processes may be related and contribute to the reduced foaling rate observed in mares with endometrosis. Our findings revealed that, despite endometrosis being a progressive condition, the myometrial expression of genes does not change progressively at different stages of this condition. However, this inconsistency in gene expression was also confirmed in a transcriptomic analysis of the endometrium at various stages of endometrosis 7 , 22 , 23 . The observed variability may be associated with several mechanisms, including the influence of the local microenvironment, cellular adaptive responses, and the activation of distinct signaling pathways and biological processes at different stages of endometrosis. Nevertheless, it is noteworthy that within the common deTARs, the direction of transcriptomic changes was consistently the same. This indicates that, while the magnitude and specific expression patterns may vary, the fundamental regulatory shifts in gene expression follow a consistent trajectory throughout the progression of endometrosis. The muscle contractions occur when myosin cross-bridges slide along actin 24 , 25 . This requires the unblocking of the myosin-binding site on actin by shifting tropomyosin with troponin 26 . In the current study, we found that the expression of TNNT1 ( troponin T1, slow skeletal type ) is up-regulated in mares with mild endometrosis. TNNT1 is annotated for i.a to muscle structure development , and regulation of muscle contraction GO BP terms. Troponin T binds the troponin components to tropomyosin 27 . The uterine troponin-tropomyosin complex inhibits uterine actomyosin ATPase 25 . A decrease in actomyosin ATPase activity can result in less effective ATP hydrolysis, leading to weaker force generation for myosin head movement and cross-bridge formation 24 . In pregnant women myometrial troponin I, skeletal, slow ( TNNI1 ) 28 , which inhibits the interaction of myosin with actin 27 , was up-regulated. Therefore, an increase in troponin level could disturb the formation of actin-myosin cross-bridges in the myometrium of mares with endometrosis. Notably, only the phosphorylated form of myosin can form cross-bridges with actin 29 . The phosphorylation status of myosin is regulated by myosin light chain kinase (MLCK) and myosin light chain phosphatase (MLCP) 29 , 30 . In the current study, we found that the expression of protein phosphatase 1 regulatory subunit 12B ( PPP1R12B ), also known as myosin phosphatase-targeting subunit 2 ( MYPT2 ) is significantly decreased in the myometrium of mares with mild endometrosis. PPP1R12B is annotated i.a. to muscle contraction , regulation of muscle contraction GO BP terms, and focal adhesion KEGG pathway and encodes both MYPT2-MBS (myosin binding subunit) which is one of the large MLCP subunits, and a small subunit of MLCP 31 , and is involved in the regulation of MLCP phosphatase activity by augmenting Ca 2+ sensitivity of the contractile apparatus 31 , 32 . For this reason, the down-regulation of PPP1R12B expression in the myometrium of mares with mild endometrosis may reduce MLCP capacity for Ca 2+ binding leading to impaired sliding of the myosin head along the actin filaments, and thus, the contractile activity of the tissue. We suggest that targeting MLCP activity might be a promising treatment for improving myometrial contractile activity in mares with endometrosis, especially in the mild stage of the condition. The intracellular Ca 2+ plays a role in a variety of molecular processes crucial for normal growth, development, and disease processes, and most importantly in the myometrium, contractile activity 33 . An elevated concentration of intracellular Ca 2+ stimulates Ca 2+ -dependent calmodulin (CAMK), activating MLCK for further myosin phosphorylation and formation of cross-bridges with actin. Therefore, impaired intracellular Ca 2+ handling may reduce the myometrial potential for generating contractile force. In the present study, we observed that genes annotated to the calcium signaling pathway , and cAMP signaling pathway in the myometrium of mares with mild and moderate endometrosis, including phospholipase C epsilon 1 ( PLCE1 ), and calcium/calmodulin-dependent protein kinase II gamma ( CAMK2G ) were down-regulated. Additionally, we determined a down-regulated expression of PLC beta 4 ( PLCB4 ) and ryanodine receptor 2 ( RYR2 ) in the myometrium of mares with moderate endometrosis (category IIB), and PLC eta 1 ( PLCH1 ) in mares with severe endometrosis (category III). PLC cleaves phosphatidylinositol 4,5-bisphosphate (PIP2) to diacylglycerol (DAG) and inositol 1,4,5-triphosphate (IP3), activating protein kinase C (PKC), and opening second-messenger-operated channels (SMOCs) in the sarcoplasmic reticulum (SR), respectively 29 , 34 , 35 . However, PLC requires the presence of Ca 2+ ions for its catalytic activity 36 . Thus, the altered expression of PLC, as well as the low basal intracellular level of Ca 2+ may interrupt further Ca 2+ release from intracellular stores. Studies on mid-pregnant rats showed that when the uterine contractions subsided, PLCB-encoding genes were down-regulated, except for PLCB4 37 . However, the basal activity of this enzyme is inhibited by ribonucleotides, including GTP-γ -S, and thus it makes it difficult to assess its activity 36 , and may also have a minor effect on uterine contractions compared to other PLCB isoenzymes 37 . Notably, this process is under the control of the RYR 33 . Studies on women and rats showed that increased expression of RYR2 is detected in the myometrium during labor and is needed for effective uterine contractions 38 – 40 . The current study showed that RYR2 in the myometrium of mares with moderate endometrosis is significantly decreased, which may contribute to the impaired contractile activity of the tissue in affected mares. Nevertheless, the lack of reproductive history of mares used in the current study has to be kept in mind when interpreting the results, as the parity of the animals may also affect contractions-related transcriptomic changes in the myometrium. Importantly, the contractile force generation in the myometrium is much greater when coupled to extracellular Ca 2+ ions entry than Ca 2+ ions release from SR 29 . The extracellular Ca 2+ entry is controlled by receptor-operated channels (ROCs) 29 , 34 , and the transient receptor protein (TRP) ion-channel family 34 . In the current study, we showed that in mares with moderate endometrosis myometrial expression of glutamate ionotropic receptor NMDA type subunit 2A ( GRIN2A ), which is part of the NMDA (N-methyl-D-aspartate) receptors assigned to ROCs 29 , 34 , was down-regulated. Also, we noted myometrial down-regulation of transient receptor potential cation channel subfamily C member 4 ( TRPC4 ) in mares with mild and severe endometrosis. Studies on women have shown that myometrial cells with depletion of TRPC4 have significantly lower extracellular calcium-dependent increases in Ca 2+ concentration. All of the above may indicate that in mares with moderate and severe condition impairment of myometrial contractility is related to the extracellular Ca 2+ entry and intracellular Ca 2+ handling. Therefore, we suggest that future therapies aiming to improve myometrial activity in mares with moderate and severe endometrosis should focus on regulating ROCs, TRP family, PLCs, and RYR2, rather than MLCK/MLCP, as valid for mild stage of the condition. Nevertheless, this requires further in-depth functional studies. Noteworthy, intracellular Ca 2+ handling and myometrial contractions are under the control of different factors, including oxytocin, prostaglandins, and cytokines 41 – 44 . This study evaluated that in the myometrium of mares with endometrosis, the expression of genes annotated to the oxytocin-signaling pathway is altered. Oxytocin was previously reported to induce myometrial contractions in reproductively healthy mares despite their age, as well as in older individuals susceptible to delayed uterine clearance 8 . Notably, Rigby et al. (2001) showed that OT-dependent increase in myometrial contractions is greater in older mares compared to younger individuals, and similar in mares susceptible to delayed uterine clearance 8 . These results suggest that alterations in contractile activity may not be strictly related to intracellular Ca 2+ concentration, but age-related changes in receptor-second messenger signaling mechanisms downstream of intracellular Ca 2+ release 8 . Among genes that were annotated to oxytocin-signaling pathway , in the myometrium of mares with different stages of endometrosis we observed increased expression of protein kinase C zeta ( PRKCZ ), and in mares with the moderate condition, protein kinase C gamma ( PRKCG ). Both genes encode different isoforms of PKC. Studies on women in late pregnancy showed that PKC zeta , by affecting actin organization and its myosin-binding properties in myometrial myocytes, stimulates myometrial contractions 45 . However, this process requires prior stimulation with endothelin-1 (EDN1) 45 , which in mares with mild endometrosis was down-regulated. Thus, a lowered expression of PLC family-coding genes, as well as RYR2 in the myometrium, may inhibit the release of Ca 2+ ions from SR, but OT may augment PLC and PKC activity for further increase of contractile force in the myometrium of mares with endometrosis. Myometrial contractions are also under the control of PGs 41 , 43 . In the current study, we showed that in mares with mild and severe endometrosis myometrial expression of prostaglandin F receptor ( PTGFR ) decreases, which may lead to decreased PGF2α-dependent release of Ca 2+ ions in myocytes. On the other hand, in mares with moderate endometrosis, we documented up-regulation of prostaglandin-endoperoxide synthase 1 ( PTGS1 ), coding a major enzyme in the prostaglandin synthesis pathway. However, the profile of myometrial prostaglandins release in mares with different severity of endometrosis was not evaluated yet, and it is not known whether stimulation with PGF2α may affect myometrial contractions of mares with endometrosis. This might be an interesting direction for further studies. The intracellular Ca 2+ concentration in the myocytes, and thus, myometrial contractions, may differ also in response to interleukin 33 (IL33) 44 , 46 . Present study showed that the expression of IL33 is lowered in the myometrium of mares with endometrosis (all stages of the condition). Consequently, it may cause a decreased influx of Ca 2+ ions to myocytes, deteriorating the contractile activity of the myometrium. Notably, IL33 is one of 200 commonly differentially expressed genes in the myometrium of mares in the course of the condition. Many of these common DEGs were assigned to adaptive immune response , antigen binding , and immunoglobulin complex GO BP terms. Importantly, immune agents including cytokines may direct myometrial fate to a contractile or quietened state 47 . Notwithstanding, all of the above-described transcriptomic changes indicate that intracellular Ca 2+ handling, including Ca 2+ influx and management from intracellular storages, may be impaired in the myometrium of mares with endometrosis, which should be evaluated in further in-depth studies. Also, we believe that evaluating the role of immune processes in these mechanisms is an interesting idea for further studies. Summarizing, our results show that the effectiveness of OT stimulation in mares with endometrosis may rely on myometrial PLC and PKC expression and activity. Also, the stimulatory effect of PG treatment on myometrial contractions might be diminished in mares with endometrosis due to decrease of PTGFR expression. Finally, we believe that future therapies for improvement of myometrial function in mares with endometrosis should focus on the involvement of immune mechanisms in the tissue, especially related to IL33 action. Based on the current study it is suggested that the observed changes in gene expression, especially in mares with endometrosis, may explain the reported histological alterations in myometrium and ECM destabilization. It is worth highlighting that in the group evaluated in the current study common DEGs with the highest fold change in the myometrium of mares with endometrosis are aggrecan ( ACAN ), assigned to the extracellular matrix, collagen-containing extracellular matrix, and extracellular matrix structural constituent GO CC terms. The expression of ACAN was lowered in the myometrium of mares with endometrosis, and based on studies in cartilage and brain, is involved in the regulation of ECM formation and tissue stiffness 48 , 49 . Notably, in the group of top-up-regulated common DEGs was classified gene encoding enzyme targeting ACAN, i.e. hyaluronidase 4 ( HYAL4 ). It is possible that both decreased expression of ACAN and increased expression of HYAL4 may contribute to the incidence of lacunae in the myometrium of affected mares, which was observed previously 13 . In mares with mild endometrosis many genes with lowered expression in the myometrium were annotated to muscle structure development , muscle tissue development , and muscle organ development GO BP terms. Among these are myocardin ( MYOCD ), and early growth response 1 ( EGR1 ). Myocardin is a coregulator of the master smooth muscle transcription factor serum response factor (SRF), regulating smooth muscle cell development and differentiation 50 , and EGR1 is predicted to regulate the second wave of gene expression inducing long-term effects on cellular growth and differentiation 51 . Based on a study in women, the MYOCD expression level in the myometrium is significantly elevated in term-pregnant females compared to non-pregnant females 52 , and it is believed that myocardin and progesterone together mediate uterine rearrangement during pregnancy 52 . As highlighted above, the reproductive history of mares used in the current study is unknown. For this reason, the observed transcriptional changes in the myometrium of mares with endometrosis might be different in maiden and multiparous animals. Nevertheless, a decreased expression of both MYOCD and EGR1 may lead to the atrophy of myocytes observed in the uterus in mares with endometrosis, even though further in-depth studies are required. The current study has shown that in the myometrium of mares with endometrosis is observed lowered expression of genes enriching ECM-receptor interaction and focal adhesion KEGG pathways , i.a. collagen IV A1 ( COL4A1 ) and integrin A1 ( ITGA1 ). Notably, we also documented in the myometrium of studied mares lowered expression of COL4A1 , COL4A5 , COL5A2 , and COL8A1 (IIA vs I), COL2A1 , COL9A1 , COL9A3 , and COL12A1 (IIB vs I), and COL28A1 (III vs I), and an increased expression of COL17A1 (IIA vs I), and COL19A1 (IIB vs I). An impairment of collagen-coding gene expression in the myometrium in mares with endometrosis may contribute to observed previous atrophy of myocytes and distension of the intercellular spaces and thus, decreased cell adhesions in this condition 13 . According to Hanada et al. 13 , an increased collagen deposition in mares with endometrosis might be connected with aging. Nevertheless, also in younger mares with developed endometrosis structural changes in the myometrium are observed 13 . The expression of collagen-coding genes changes in the course of pregnancy and is crucial for switching myometrial phenotype from contractile to quietened one 53 . Notably, we showed increased myometrial expression of genes coding for matrix metallopeptidases (MMPs), specifically MMP9 in IIA vs I, and MMP15 in III vs I. A study on women shows that uteroplacental expression of MMP9 is lowered in physiological vs preeclamptic pregnancy 54 . Furthermore, a study on rats shows that in pregnant vs virgin females uterine expression of MMP9 is increased, and inhibition of its activity increases myometrial contractions 55 . Thus, in mares with endometrosis, similar mechanisms may occur as such, during aging or pregnancy, leading to structural changes and quietened myometrial contractility. Providing that in each stage of endometrosis progression we observed alterations in the myometrial expression of genes encoding ECM components as well as disturbance in ECM-receptor interactions and focal adhesion, we are positive that this is the future of designing therapeutic strategies for mares with this condition. ECM-related changes may impact myometrial cells fate and direct switching phenotype to quietened one, impairing function of the tissue. The occurrence of such changes may further aggravate the condition. Thus, further in-depth functional studies are required for effective manipulation of ECM composition in the myometrium of mares with endometrosis.

The transcriptomic activity of the myometrium is altered in the progression of a mare’s endometrosis. The most abundant changes were observed in the myometrium of mares with mild endometrosis compared to moderate and severe stages of this condition. Evaluated changes in gene expression may indicate impaired function of contractile machinery, mechanisms regulating Ca 2+ influx and handling, and changes in ECM composition in myometrium. All of the above may contribute to observed previously histological alterations and deteriorated contractile activity of the tissue in mares with endometrosis. Based on the provided results, changes occurring in the myometrium may contribute to the development of persistent endometritis, delayed uterine clearance, and consequently endometrosis but they may also be a consequence of the pathological processes occurring in the endometrium. This is yet to be evaluated. However, since alterations in the myometrium are likely to contribute to the impairment of reproductive potential in affected mares, further studies must be carried out to unravel the intricate cross-talk between mare endometrium and myometrium.

Endometrosis is a chronic degenerative uterine condition characterized primarily in horses ( equine ) that, among others, is distinguished by fibrosis around the endometrial glands and in the endometrial stroma. This term was introduced by Kenney 1 to define changes in the mare uterus previously referred as chronic degenerative endometritis. In horses, as the severity of endometrosis increases, the expected foaling rate decreases by 20–50%, 50–90%, and more than 90% for mares with endometrium category IIA (mild condition), IIB (moderate condition), and III (severe condition), respectively, in comparison to mares without endometrial alterations (category I) 2 . Thus, endometrosis substantially reduces the reproductive potential of mares causing a considerable economic problem for horse breeders 2 . Endometrosis is not only limited to the histological changes in the endometrium but also encompasses changes in the endometrial microenvironment 3 – 6 . These changes include alterations in the endometrial expression of steroid hormone receptors 3 , and prostaglandin (PG) synthases 6 , as well as endometrial secretion of i.a. uteroglobin, uterocalin, calbindin, and glycogen 3 . Furthermore, it has recently been determined that the endometrial transcriptome profile varies in mares with different stages of endometrosis 7 . Therefore, extensive research is currently being conducted to enhance our comprehension of the etiology and pathogenesis of equine endometrosis to develop effective treatments for this condition and prevent impaired reproductive potential in mares. Endometrosis may develop from persistent endometritis caused by delayed uterine clearance 8 , 9 . This may result from impaired function of the myometrium, however, this was poorly investigated so far, especially on the molecular level. It is clear that the myometrium’s function is vital for the transportation of sperm and embryos, as well as delivery 10 – 13 , thereby playing a significant role in preserving the female reproductive capacity. Moreover, in addition to the contractile activity, mare’s myometrium in an endocrine organ producing many regulatory factors, including PGs 14 . The atrophy of the smooth muscle fascicles and cells in the uterine muscle layer, as well as fatty degeneration of atrophic myocytes and an increase in the deposition of collagen and elastin fibers between the smooth muscle bundles have been documented in cases of endometrosis in mares 15 . Moreover, it is well known that in mares susceptible to persistent endometritis a lowered contractile activity of the myometrium in response to inflammation is observed 16 . In consequence, it may lead to a deteriorated uterine clearance and the development of endometrosis 9 . Nevertheless, as reported by Rigby et al. 8 , the myometrial contractions may be increased in response to oxytocin (OT) and PGs even in older mares susceptible to delayed uterine clearance, but the mechanism of this physiological effect is yet to be cleared 8 . Early works have proposed, as a possible cause of impaired uterine clearance, the loss of structural support of the caudal reproductive tract and stretching of the broad ligaments in old multiparous mares, which may cause uterus dropping and tilting ventrally in relation to the pelvic brim 17 , 18 . The molecular mechanisms responsible for the impaired myometrial function in susceptible mares with persistent endometritis and/or endometrosis are complex and poorly understood. Therefore, we aimed to determine the transcriptional profile of the myometrium of mares with different endometrial scores (categories IIA, IIB, and III according to the Kenney and Doig classification system 2 ) compared to mares without changes in the endometrium (category I). RNA-seq was used to identify differentially expressed transcript active regions (deTARs) which correspond to protein-coding transcripts, i.e. differentially expressed genes (DEGs) and non-protein coding, uncovered regulatory region groups, in the myometrium of mares with endometrosis, in contrast to endometria without changes, and to investigate potential changes in biological processes and signaling pathways that may affect mares’ uterine functions. To the best of our knowledge, this is the first study that provides a transcriptomic characterization of the myometrium of mares during the mid-luteal phase of the estrous cycle in the progression of endometrosis. Transcriptomic analysis can help fill the gaps in our understanding of myometrial changes and identify novel pathways and regulators present in the myometrium at different stages of endometrosis.

Supplementary Information 1. Supplementary Information 2. Supplementary Information 3. Supplementary Information 4. Supplementary Information 1. Supplementary Information 2. Supplementary Information 3. Supplementary Information 4.