Transcriptomics of receptive endometrium in women with sonographic features of adenomyosis

Background: Women with uterine adenomyosis seeking assisted reproduction have been associated with compro- mised endometrial receptivity to embryo implantation. To understand the mechanisms involved in this process, we aimed to compare endometrial transcriptome profiles during the window of implantation (WOI) between women with and without adenomyosis.

We obtained endometrial biopsies LH-timed to the WOI from women with sonographic features of adenomyosis (n=10) and controls (n=10). Isolated RNA samples were subjected to RNA sequencing (RNA-seq) by the Illumina NovaSeq 6000 platform and endometrial receptivity classification with a molecular tool for menstrual cycle phase dating (beREADY®, CCHT). The program language R and Bioconductor packages were applied to analyse RNA-seq data in the setting of the result of accurate endometrial dating. To suggest robust candidate pathways, the identified differentially expressed genes (DEGs) associated with the adenomyosis group in the receptive phase were further integrated with 151, 173 and 42 extracted genes from published studies that were related to endometrial receptivity in healthy uterus, endometriosis and adenomyosis, respectively. Enrichment analyses were performed using Cytoscape ClueGO and CluePedia apps.

Out of 20 endometrial samples, 2 were dated to the early receptive phase, 13 to the receptive phase and 5 to the late receptive phase. Comparison of the transcriptomics data from all 20 samples provided 909 DEGs (p<0.05; nonsignificant after adjusted p value) in the adenomyosis group but only 4 enriched pathways (Bonferroni p value < 0.05). The analysis of 13 samples only dated to the receptive phase provided suggestive 382 DEGs (p<0.05; non- significant after adjusted p value) in the adenomyosis group, leading to 33 enriched pathways (Bonferroni p value < 0.05). These included pathways were already associated with endometrial biology, such as “Expression of interferon (IFN)-induced genes” and “Response to IFN-alpha” . Data integration revealed pathways indicating a unique effect of adenomyosis on endometrial molecular organization (e.g., “Expression of IFN-induced genes”) and its interfer- ence with endometrial receptivity establishment (e.g., “Extracellular matrix organization” and “Tumour necrosis factor production”).

Accurate endometrial dating and RNA-seq analysis resulted in the identification of altered response to IFN signalling as the most promising candidate of impaired uterine receptivity in adenomyosis. © The Author(s) 2021. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. The Creative Commons Public Domain Dedication waiver (http:// creat iveco mmons. org/ publi cdoma in/ zero/1. 0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Open Access *Correspondence: [email protected] 5 Department of Gynaecology, University Medical Centre Maribor, 2000 Maribor, Slovenia Full list of author information is available at the end of the article Page 2 of 16Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 Executive summary of the study 1. Adenomyosis has been associated with lower preg - nancy rate in infertility treatments. 2. Molecular knowledge of endometrial receptivity in women with adenomyosis is scarce and limited to studies with selected candidate genes and one genome-wide study performed using microarrays. 3. Therefore, we performed the first transcriptome sequencing of endometrial samples LH-timed to the expected WOI (LH+7 – LH+9) between women with (n = 10) and without (n = 10) sono - graphic features of adenomyosis. 4. Transcriptomics data comparison of 10 adeno - myosis cases and 10 normal controls provided 909 DEGs (p<0.05; nonsignificant after adjusted p value), but functional enrichment analysis iden - tified 4 pathways (Bonferroni p value < 0.05) that were not directly associated with endometrial biol - ogy. 5. Retrieved endometrial biopsies were applied for the external molecular tool beREADY ® (CCHT, Estonia) to verify their receptivity status on the basis of the gene expression signature associated with endometrial receptivity. Out of 20 samples, 2 were classified as early receptive, 13 as receptive and 5 as late receptive. 6. Two early- and 5 late-receptive samples were excluded from the RNA-seq dataset to prevent the impact of early- and late-secretory phases of the menstrual cycle on transcriptomics analysis associated with endometrial receptivity. The RNA- seq dataset of the remaining 8 adenomyosis cases and 5 control receptive endometrial samples was reanalysed, and 382 DEGs (p < 0.05; nonsignificant after adjusted p value) were identified, resulting in 33 enriched pathways (Bonferroni p value < 0.05) that have already been associated with endometrial biology. 7. The 382 identified DEGs were further integrated with the most extensive set of genes from the lit - erature associated with endometrial receptivity in the healthy uterus, endometriosis (model disease to study persistence of gynaecological pathology on endometrial molecular organization) and adeno - myosis to provide candidate pathways characteriz - ing the role of adenomyosis on endometrial molec - ular organization. 8. Integrative enrichment analysis provided candi - date pathways that may indicate a unique effect of adenomyosis on endometrial molecular organiza - tion (e.g., “Expression of IFN-induced genes”) and its interference with endometrial receptivity estab - lishment (e.g., “Extracellular matrix organization” , “Tumour necrosis factor production” and “Regula- tion of reproductive process”). 9. Identification of robust endometrial pathways and associated genes could lead to the development of molecular tools for endometrial receptivity exami - nation that would be specific for women with adenomyosis. 10. Accurate endometrial receptivity examination in infertile adenomyosis patients could better ver - ify whether endometrial-associated factors are a source of recurrent implantation failures.

Adenomyosis is a common acquired uterine anomaly characterized by the presence of endometrial glands and stroma within the myometrium. Advances in imag - ing techniques in the last decade have enabled the diag - nosis of adenomyosis [1] in a large proportion of women undergoing infertility diagnostics [2, 3]. Since subtle sonographic signs of adenomyosis are becoming easier to recognize, adenomyosis is diagnosed with increasing frequency. Previous retrospective studies have shown the association between adenomyosis and lower embryo implantation rates and higher miscarriage rates [4–6]. Several functional and molecular aberrations could be responsible for altered endometrial receptivity to embryo implantation and lower fecundity in women with adeno - myosis. It has been suggested that the disruption of the junctional zone architecture by adenomyosis could lead to altered contractility and interrupt endometrial recep - tivity [7,  8]. Other suggested causes affecting endome - trial receptivity in women with adenomyosis could be increased levels of oxidative stress [9–11], abnormal endometrial vascularity [12, 13] and functional disorgan- ization at the molecular level [14–17]. In our previous study [18], we gathered proteins, genes and functional noncoding RNAs (ncRNAs) shown to be dysregulated in the endometrium of women with adeno - myosis during the expected window of implantation (WOI). Bioinformatics approaches were used to integrate

Adenomyosis, Assisted reproductive techniques (ART), Data integration, Endometrial receptivity, Enrichment pathway analysis, Omics approaches, RNA-seq, Systems biology, Transcriptomics, Window of implantation Page 3 of 16 Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 retrieved loci with endometrial receptivity genes from the literature associated with healthy (normal) uteri to iden - tify candidate dysregulated mechanisms involved in the regulation of embryo implantation in adenomyosis. In addition, we included better characterized endometriosis as a model disorder to study the impact of gynaecological pathology on endometrial molecular organization [18]. Numerous published genome-wide studies associated with the endometrial molecular background in women with endometriosis enabled us to develop a catalogue of genes sorted according to the phases of the menstrual cycle [19]. Genes sorted in the mid-secretory phase corre- sponding to the appearance of the WOI were used for the integrative analysis mentioned above [18]. The identified enriched “Signalling by interleukins” and “Interleukin-4 and interleukin-13 signalling” pathways were prioritized, and the corresponding mapped LIF, SOCS3, IL10, IL6, JUNB and FOS genes were validated. Since downregu - lated expression levels of selected genes in adenomyosis compared to the control group showed no statistical sig - nificance, we assumed that comprehensive endometrial transcriptomics profiling would be an appropriate next step to identify adenomyosis-specific loci [18]. To date, there is only one transcriptomics study [20] profiling the endometrium in the expected WOI using microarrays, which identified 34 differentially expressed genes (DEGs) in women with adenomyosis wishing to conceive compared to healthy women [20]. The meth - odological improvement of transcriptome profiling from hybridization-based microarrays to next-generation sequencing (NGS) platforms provides more comprehen - sive insight into expression signatures and enables iden - tification of minor differences between study groups [21]. Millions of reads generated by RNA sequencing (RNA- seq) can be aligned to a reference genome, reference transcripts or references assembled de novo for the entire transcriptome to be surveyed. Thus, additional biological constituents can be identified, and a more precise assess - ment of transcript expression levels can be obtained [22]. The first aim of this study was to perform RNA-seq of endometrial samples dated to the WOI between women with and without sonographic features of adenomyosis to identify DEGs. The second aim was to perform enrich - ment analysis of identified DEGs alone and together with endometrial receptivity genes from the literature to pro - vide robust candidate pathways related to altered molecu- lar background of endometrial receptivity in adenomyosis.

Study cohorts We designed a prospective observational study includ - ing women scheduled for medically assisted reproduc - tion at the Department of Reproductive Medicine and Gynaecological Endocrinology, University Medical Centre Maribor, Slovenia between 2018 and 2020. The inclusion criteria were as follows: age ≤ 42 years, regular menstrual cycle 24 – 36 days in length, no cur - rent hormonal treatment, controlled ovarian stimula - tion (COS), ovulation triggering or vaginal progesterone for luteal support at least two months prior to endo - metrial biopsy. The exclusion criteria were anovulatory menstrual cycles, polycystic ovary syndrome (PCOS), previous surgical treatment of endometriosis or uterine surgical procedures, sonographic evidence of fibroids, endometrial polyps, hydrosalpinges, and evidence of ovarian or deep infiltrating endometriosis (unless other - wise noted in Table 1). In our clinic, all women undergo - ing assisted reproductive techniques (ART) have a prior transvaginal ultrasound (TVUS) examination, typically performed in the proliferative phase of the menstrual cycle. Women with echographic evidence of adenomyo - sis were considered eligible for the study, and the con - trol group was composed of women with normal uteri seeking ART due to male or tubal factors of infertility. On the day of endometrial sampling, all women under - went TVUS performed by a single expert sonographer (level 3 according to European Federation of Societies for Ultrasound in Medicine and Biology). In all women, comprehensive 2-D and 3-D ultrasound using high-range equipment was performed with a 10  MHz transvaginal transducer (Voluson E8 Expert, GE Health care, Austria GmbH & Co OG, Zipf, Austria). Diagnostic criteria for adenomyosis were based on previously published criteria [23]. The diagnosis of adenomyosis was confirmed when one of the following sonographic criteria was met: asym - metrical myometrial thickening not caused by the pres - ence of fibroids, linear endometrial striations, irregular endometrial-myometrial junction, parallel shadowing, or the presence of myometrial cysts or hyperechoic islands [23]. Adenomyosis was classified as mild by subjective assessment, but in general, it was assessed in line with previously described principles. This was when only focal areas of adenomyosis were seen or when adenomyosis was present only in the inner third of the myometrium [24]. Demographic and clinical characteristics of partici - pants, including age, body mass index (BMI), endometrial thickness at the time of endometrial biopsy and the num- ber of previous ART cycles, are presented as the median (range) and were compared between study groups using the nonparametric Mann–Whitney U-test in SPSS 25.0 software (IBM Corporation, Armonk, NY, USA). Statisti - cal significance was set at p value < 0.05. Endometrial sample collection Endometrial biopsy sampling was conducted in a natu - ral menstrual cycle, and women were scheduled for Page 4 of 16Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 cycle monitoring by urinary luteinizing hormone (LH) tests (Hangzhou AllTest Biotech Co., Ltd, Hangzhou, P .R. China). Women were scheduled for endometrial sampling conducted by the Pipelle endometrial suc - tion curette (the Probet, Gynetics Medical Products N.V., Lommel, Belgium) in the expected WOI on the day between LH+7 to LH+9 after a participant’s LH surge determination (day LH+0). Retrieved endometrial samples were immediately placed in RNAlater solution (Thermo Fisher Scientific Baltics UAB, Vilnius, Lithu - ania), stored overnight at +4 °C and then transferred to –80 °C until RNA isolation was performed. Total RNA isolation and quality control Total RNA was isolated using the miRNeasy Mini Kit (Qiagen GmbH, Hilden, Germany) according to the man- ufacturer’s instructions. Each whole-tissue endometrial sample was first disrupted with a Bullet Blender Storm Pro homogenizer (Next Advance, lnc., Troy, NY, USA) using 1 mm zirconium oxide beads in 700 µL of QIAzol Lysis Reagent from the miRNeasy Mini Kit. After 5 min of incubation at room temperature, 140 µL of chloro - form was added to the homogenate, and the solution was shaken vigorously. The sample was then centrifuged at 12 000 rfc for 15 min at 4 °C. The upper aqueous phase (approximately 300 µL) was transferred to a new Eppen - dorf tube, and 1.5 volumes of ethanol were added. The samples were then pipetted to RNA binding miRNeasy Mini spin columns and washed using RWT Buffer and RPE Buffer solutions of the miRNeasy Mini Kit. Total RNA was eluted in 50 µL of RNase-free H2O. The quantity and purity of each RNA sample were assessed with Synergy 2 spectrophotometric measure - ments (BioTek Instruments, Winooski, VT, USA). RNA integrity number (RIN) was estimated on the 2100 Bio - analyser system (Agilent Technologies, Waldbronn, Germany) using the RNA Nano 6000 Assay Kit (Agilent Technologies, Waldbronn, Germany). After passing those quality controls, each RNA sample was used for cDNA library construction and subsequent RNA-seq and for accurate endometrial dating of retrieved biopsies. Accurate endometrial dating One part of each RNA sample was shipped on dry ice to the Competence Centre on Health Technologies, CCHT, Tartu, Estonia, where endometrial receptiv - ity testing was performed using the beREADY ® test [25] (https:// berea dy. ccht. ee/). Endometrial dating was performed according to the established protocol using targeted allele counting by sequencing (TAC-seq) meth - odology [26] to explore the expression levels of 57 well- described endometrial receptivity genes [27]. The results of the beREADY ® test were provided in five phases: “pre-receptive” , “early-receptive“, “receptive” , “late-recep- tive” , and “post-receptive” . The purpose of endometrial dating was to accurately classify the receptivity status of LH-timed biopsies to remove samples that could lead to possible biases in gene expression analysis associated with endometrial receptivity in adenomyosis. Library preparation and RNA‑seq Both lncRNA and mRNA 150 bp paired-end libraries were constructed and subsequently sequenced by Novogene Bioinformatics Technology Co., Ltd. (Hong Kong, China). Briefly, a total amount of 2 µg of RNA per sample was used for cDNA sequencing library preparation. Ribosomal RNA (rRNA) was removed using the Epicentre Ribo-zeroTM rRNA Removal Kit (Epicentre, Brooklyn, NY, USA), and the remaining RNA was used for library generation by the NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA). First, rRNA-depleted RNA samples were fragmented followed by first- and second-strand cDNA synthesis. The sequencing adaptors were ligated, and library fragments were purified to obtain cDNA fragments 150~200 bp in length. Polymerase chain reaction (PCR) amplification of size-selected, adaptor- ligated cDNA was performed using universal PCR primers and index primers. Index-coded samples were clustered by Illumina TruSeq PE Cluster Kit v3-cBot-Hs. Libraries were sequenced on an Illumina NovaSeq 6000 platform, which generated 150 bp paired-end reads. RNA‑seq data alignment and identification of DEGs Raw sequence reads were trimmed by Novogene in- house Perlscript to remove raw reads with adapter con - tamination and reads containing poly-N and low-quality reads. The RNA-seq data presented in this study are deposited in the Gene Expression Omnibus (GEO) data - base with accession number  GSE185392. Provided raw fastq files were first evaluated with FastQC v.0.11.9 soft - ware (http:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/) to obtain a quality profile of the reads. The statistical environment R v.4.0.2 (R Core Team 2020, Vienna, Austria) and contributed packages from the R software repository Bioconductor (http:// www. bioco nduct or. org/) were used for high-throughput sequence data analysis. Raw paired-end reads were aligned to the UCSC Homo sapiens hg19 reference genome using the Rsubread v.2.2.4 R package [28, 29]. Properly mapped reads were sorted in files with binary alignment/map (BAM) format. Mapped reads were counted and assigned to genomic features using fea - tureCounts [30] with the requirement that both ends should be mapped. Counts per million (CPMs) were calculated using the edgeR v.3.30.3 R package [31]. Genes expressed at low levels were filtered out based on Page 5 of 16 Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 CPMs corresponding to read counts of 10, and retained genes were normalized using the trimmed mean of M values method (TMM) [32]. Subsequently, mean-vari - ance modelling at the observational level transforma - tion (VOOM) was applied [33]. Differential expression analysis of the adenomyosis group relative to the con - trol group was determined in two RNA-seq datasets using linear models and empirical Bayes implemented in the limma v.3.44.3 R package [34]. RNA-seq datasets were composed of libraries on the basis of the results of endometrial dating of corresponding samples. The first dataset contained all LH-timed samples, while the sec - ond dataset contained only samples dated to the recep - tive phase. Differential expression was considered for genes with a p value < 0.05 regardless of the adjusted p value obtained after multiple testing corrections. Integration of identified DEGs in the adenomyosis group with endometrial receptivity genes from the literature Identified DEGs between adenomyosis cases and con - trols using samples dated to the receptive phase were applied for integrative bioinformatics analysis to pro - vide robust candidate pathways associated with altered molecular background of endometrial receptivity in adenomyosis. DEGs were applied for enrichment rea - nalysis with lists of 42, 173 and 151 genes associated with endometrial receptivity in adenomyosis, endo - metriosis and healthy uterus, respectively, that were retrieved from the literature in our previous study [18]. Genes associated with endometriosis presented a model to study the impact of gynaecological pathology on endometrial molecular organization. Genes asso - ciated with a healthy uterus were used as a reference molecular background required for endometrial recep - tivity establishment. Two enrichment analyses were performed using two different gene lists associated with adenomyosis. The first adenomyosis gene list con - tained only 382 DEGs of the present sequencing experi - ment, while the second list combined 382 DEGs with 42 genes from the literature (in total, 424 genes). The first enrichment analysis was performed by integrating the adenomyosis gene list with 382 DEGs, the endometrio - sis list with 173 genes and the healthy uterus list with 151 genes. Second, enrichment analysis was performed using adenomyosis, healthy uterus and endometriosis lists with 424, 151 and 173 genes, respectively. The gene lists used are provided in Additional file 1 . Functional enrichment analyses DEGs (p<0.05) that were identified by transcriptom - ics data comparison of endometrial samples between adenomyosis cases and controls were subjected to functional enrichment analyses using ClueGO v.2.5.8 [35] and CluePedia 1.5.8 [36] apps of Cytoscape v.3.8.2 software [37]. The same bioinformatics tools were used for enrichment analyses employing integrated gene lists associated with adenomyosis, endometriosis and healthy uterus. When analysing identified sets of DEGs associated with the present adenomyosis groups, up- and down - regulated genes were separately uploaded as two clus - ters in the ClueGO app, which gave a unique colour marker to each gene set. When performing enrich - ment analyses of integrated gene lists associated with different gynaecological conditions, each gene list was uploaded as a cluster in the ClueGO app to distinguish study groups according to colour markers of the cluster. Each enrichment analysis was applied by repre - sentative Gene Ontology Biological Process (GO_BP), Reactome Pathways and Reactome Reactions ontolo - gies. Only enriched pathways (Reactome pathways/ reactions and GO_BP terms) with corrected p values < 0.05 according to the Bonferroni step down test were considered. The identified pathways were sorted into groups based on their common biological role and associated genes (kappa score) and further projected into functionally organized networks. The size of nodes in the generated networks was correlated with the obtained p value. The pathway with the highest signifi - cant value was considered to be the leading term of a group and was therefore highlighted in the network by a large name label and a statistical summary. The Clue - Pedia app was further applied to visualize shared initial genes within or between functional network groups. The proportion of visible genes mapped to each path - way was also determined. When more than 60% of mapped genes originated from one of the clusters, a pathway was shown in the network with the predefined colour of this cluster.

An overview of the study is outlined in Fig. 1 . Participant characteristics The demographic and clinical characteristics of the study cohorts are summarized in Table 1 . Total RNA quality Total RNA was isolated from 20 endometrial samples, 10 from the adenomyosis group and 10 from the con - trol group. The A260/A280 ratios and RIN values of all RNA samples were above 2.0 and >8.5, respectively, and were further used for endometrial dating and RNA-seq. Page 6 of 16Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 Endometrial receptivity classification of LH‑timed biopsies The results of endometrial receptivity testing performed on each endometrial RNA sample are provided in Table 2. According to the test, 13 out of 20 samples were classified in the receptive phase (8 adenomyosis cases and 5 controls), 2 samples in the early receptive phase, 5 samples in the late receptive phase and zero samples in the pre- or postreceptive phases. The timing of endometrial biopsy and measured receptivity status are provided for each sample, followed by a summary of mapped RNA-seq reads and library size after filtering for low gene expression. Abbrevia - tions “ A” refer to adenomyosis and “K” to control sam - ples. Timing of biopsy refers to the day after luteinizing hormone (LH) peak determination (LH+ 0) by urinary LH test. RNA‑seq outcome parameters RNA-seq of endometrial samples utilizing the Illumina NovaSeq 6000 platform generated between 53,026,608 and Fig. 1 The overview of the present study. Only ovulatory women with regular menstrual cycles were included. Endometrial biopsies, 10 in the adenomyosis group and 10 the in control group, were conducted between 7 and 9 days post urinary LH peak corresponding to the expected WOI. Isolated RNA samples were used for RNA-seq using the Illumina NovaSeq 6000 platform and for accurate endometrial receptivity classification using the beREADY® molecular test. When analysing transcriptomics data of all 20 LH-timed biopsies, 909 DEGs (p < 0.05, nonsignificant after corrected p value) were associated with the adenomyosis group. Downstream functional enrichment analysis of these genes identified no strong candidate mechanisms associated with endometrial molecular biology. According to endometrial receptivity testing, 2 out of 20 samples were classified as early-receptive, 13 as receptive and 5 as late-receptive during the menstrual cycle. To prevent early- and late-secretory phases of the menstrual cycle on the transcriptomics analysis associated with endometrial receptivity, samples dated to the early- and late receptive phases were omitted from the RNA-seq dataset. The remaining transcriptomics data of 8 adenomyosis cases and 5 control samples dated to the receptive phase were reanalysed. The 382 identified DEGs (p <0.05, nonsignificant after corrected p value) in the adenomyosis group were further enriched in more robust candidate pathways, including “Expression of IFN-induced genes” , “Response to interferon-alpha” and “ISG15-protein conjugation” . The 382 identified DEGs were further integrated with 42, 173 and 151 genes from the literature associated with endometrial receptivity in adenomyosis, endometriosis and healthy uterus, respectively, to propose a molecular background of endometrial receptivity under adenomyosis. Abbreviations: DEGs = differentially expressed genes; LH = luteinizing hormone; RIN = RNA integrity number; WOI = window of implantation. The images of the NovaSeq 6000 Sequencing System and Gynetics suction curette were obtained from official pages Illumina.com and gyneatic.com, respectively Page 7 of 16 Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 75,561,205 reads per sample, with an average of 64,266,576 reads. Quality analysis of raw RNA-seq reads by FastQC revealed that each fastq file contained reads 150 base pairs (bp) in length with a mean per base sequence quality score (Phred score) of 36 and thus each file was considered for downstream bioinformatics analysis. Table 2 summarizes the number and proportion of mapped raw reads to the hg9 reference genome for each sample and obtained library sizes after filtering low-expression genes. Identified DEGs associated with adenomyosis group Differential expression analyses were conducted using two RNA-seq datasets constructed of samples according Table 1 Adenomyosis and control group characteristics The median (range) is indicated for age, BMI, endometrial thickness and number of performed ART cycles (in vitro fertilization (IVF) and/or intracellular sperm injection (ICSI) treatments). P values are based on Mann–Whitney U test. Primary sterility refers to women who have never been pregnant (nulligravid), and secondary sterility refers to women who have already achieved pregnancy (gravida) or delivery (parous). Abbreviations: BMI = body mass index; ART = assisted reproductive technique Characteristic Adenomyosis group (N = 10) Control group (N = 10) p value Age (years) 35 (30–39) 34.5 (30–42) 0.621 BMI (kg/m2) 27.2 (17.8–34.6) 21 (17.3–30.1) 0.112 Endometrial thickness (mm) 7.1 (4.6–11.2) 8 (6.2–10.1) 0.082 Number of performed ART cycles 2 (1–4) 4 (1–6) 0.353 Women sterility status: Primary sterility (nulligravid) 6 5 Secondary sterility (gravida or parous) 4 5 Factor of infertility: Male 5 8 Tubal 2 1 History of endometriosis 1 0 Idiopathic infertility 2 1 Table 2 Characteristics of endometrial RNA samples used in the study Sample ID Day of biopsy sampling Endometrial dating by the beREADY® test Number of mapped RNA‑seq reads Proportion of mapped RNA‑seq reads Library size after normalization A10 LH+8 receptive 70,040,785 97.97% 25,194,593 A12 LH+7 receptive 71,578,907 98.25% 24,095,230 A18 LH+8 receptive 59,462,597 97.85% 19,682,295 A20 LH+7 receptive 73,911,463 97.82% 26,033,066 A21 LH+8 late-receptive 66,603,640 98.34% 24,868,485 A29 LH+7 receptive 52,378,707 98.78% 19,079,266 A31 LH+7 receptive 71,617,925 98.78% 23,269,754 A3 LH+7 receptive 60,702,730 98.69% 19,575,933 A5 LH+7 receptive 63,795,289 98.02% 21,175,249 A9 LH+7 early-receptive 64,755,616 98.12% 26,091,668 K11 LH+9 receptive 60,841,108 98.02% 22,021,130 K15 LH+8 receptive 64,147,880 98.25% 23,535,991 K17 LH+8 late-receptive 64,945,168 97.84% 21,964,495 K22 LH+9 late-receptive 66,813,165 97.91% 22,118,116 K23 LH+8 late-receptive 50,492,938 98.76% 16,086,563 K24 LH+7 receptive 58,204,527 98.72% 17,532,809 K26 LH+9 late-receptive 55,574,782 98.87% 19,198,159 K27 LH+9 receptive 64,512,592 98.81% 19,919,022 K28 LH+7 early-receptive 57,037,427 98.63% 18,335,825 K8 LH+7 receptive 66,124,560 98.08% 25360,931 Page 8 of 16Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 to the results of endometrial receptivity testing. The first RNA-seq dataset was composed of all 20 samples: 13 receptive, 2 early- and 5 late-receptive samples. The second RNA-seq dataset was composed of 13 receptive samples only, while 2 early- and 5 late-receptive samples were omitted to exclude the influence of early- and late- secretory phases of the menstrual cycle on endometrial transcriptomic analysis associated with endometrial receptivity. Transciptomics data comparison of 10 adenomyosis and 10 control samples resulted in 909 DEGs (p<0.05) associated with the adenomyosis group (the entire list of 909 DEGs is presented in Additional file  2). Accord - ing to the HUGO Gene Nomenclature Committee (HGNC) (version updated March 23, 2021) nomencla - ture system (https:// www. genen ames. org/), different locus types were identified, including 829 protein-cod - ing genes (mRNAs), 27 long noncoding RNAs (lncR - NAs), 5 microRNAs (miRNAs), 5 small nucleolar RNAs, 28 pseudogenes, 1 complex locus constituent and 14 loci that were not mapped in the HGNC data - base. Among 909 DEGs, 487 genes (452 mRNAs, 11 lncRNAs and remaining other loci types) were upreg - ulated, and 422 genes (376 mRNAs, 16 lncRNAs and remaining other loci types) were downregulated. How - ever, the fold change (FC) of expression levels between study groups was nonsignificant after the application of multiple comparison correction. Transcriptomics data comparison of 8 adenomyosis cases and 5 control endometrial samples with confirmed receptive phase provided 382 DEGs (p < 0.05) associ - ated with the adenomyosis group (the entire list of 382 DEGs is presented in Additional file  3). According to the HGNC nomenclature system, 323 loci were mRNAs, 23 lncRNAs, 21 pseudogenes, 4 miRNAs, 1 complex locus constituent, 1 T cell receptor gene and 9 uncharacter - ized. Among 382 DEGs, there were 166 upregulated (137 mRNAs, 14 lncRNAs and remaining other loci types) and 216 downregulated (186 mRNAs, 9 lncRNAs and remain- ing other loci types) genes. However, there were no sig - nificant DEGs between the study groups according to the adjusted p value. Among 382 DEGs in the adenomyosis group, up to the top 10 up- and downregulated mRNAs and lncRNAs with the highest logFC values of expression levels are presented in Tables 3 and 4, respectively. Enriched pathways associated with identified DEGs Functional enrichment analysis of 909 DEGs associated with the adenomyosis group that were obtained from transcriptomics comparison of all 20 endometrial sam - ples provided only 4 enriched GO terms sorted within 2 functionally organized network groups: “Intracellular Table 3 Top upregulated mRNAs and lncRNAs DEGs were insignificant after multiple testing correction of the p value. Abbreviation “FC” refers to fold change of expression levels ENTREZ ID HGNC symbol Long name Locus type logFC p value 259289 TAS2R43 taste 2 receptor member 43 mRNA 0.9484 0.0215 4250 SCGB2A2 secretoglobin family 2 A member 2 mRNA 0.9244 0.0150 1747 DLX3 distal-less homeobox 3 mRNA 0.9212 0.0247 54959 ODAM odontogenic, ameloblast associated mRNA 0.9184 0.0460 353091 RAET1G retinoic acid early transcript 1G mRNA 0.8803 0.0346 84072 HORMAD1 HORMA domain containing 1 mRNA 0.8581 0.0031 563 AZGP1 alpha-2-glycoprotein 1, zinc-binding mRNA 0.8259 0.0463 100507436 MICA MHC class I polypeptide-related sequence A mRNA 0.8127 0.0203 7348 UPK1B uroplakin 1B mRNA 0.8086 0.0224 158131 OR1Q1 olfactory receptor family 1 subfamily Q member 1 mRNA 0.8077 0.0466 100505967 LINC00645 long intergenic non-protein coding RNA 645 lncRNA 2.9900 0.0056 100130231 LINC00861 long intergenic non-protein coding RNA 861 lncRNA 2.0003 0.0177 654412 FAM138B family with sequence similarity 138 member B lncRNA 1.8938 0.0452 100505921 GLCCI1-DT GLCCI1 divergent transcript lncRNA 1.8543 0.0323 100506334 LINC00649 long intergenic non-protein coding RNA 649 lncRNA 1.7910 0.0314 284578 MFSD4A-AS1 MFSD4A antisense RNA 1 lncRNA 1.6098 0.0022 283876 LINC00921 long intergenic non-protein coding RNA 921 lncRNA 1.3284 0.0443 100505625 LINC02102 long intergenic non-protein coding RNA 2102 lncRNA 1.3182 0.0300 100507398 INTS6-AS1 INTS6 antisense RNA 1 lncRNA 1.2787 0.0023 93653 ST7-AS1 ST7 antisense RNA 1 lncRNA 1.2489 0.0120 Page 9 of 16 Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 lipid transport” (11 mapped genes, corrected p value 2.15 × 10-5) and “Icosanoid receptor activity” (6 mapped genes, corrected p value 2.02 × 10-5). Sorted pathways with associated genes in networks are presented in Fig.  2a. The results of the enrichment analysis are sum - marized in Additional file 4. Functional enrichment analysis of 382 DEGs associ - ated with the adenomyosis group that were obtained by transcriptomics analysis of endometrial samples in the receptive phase resulted in 33 enriched pathways, includ- ing 20 GO_BP terms, 6 Reactome pathways and 7 Reac - tome reactions. They were sorted into 7 network groups to remove redundancy, which is visualized in Fig.  2b. The highest proportion of enriched pathways was related to mechanisms of response to interferon (IFN) signal - ling, in particular antiviral response (presented in higher resolution in Fig.  2c). Most of the downregulated genes were mapped in the following network groups: “Expres - sion of IFN-induced genes” (BST2, IFI35, IFIT1, IFITM1, ISG15, MX1, OAS2, OAS3 and STAT1 were down- and IRF6 was upregulated, corrected p value 2.08 × 10-6), “Response to interferon-alpha” (BST2, EIF2AK2, IFITM1, and LAMP3, corrected p value 8.75 × 10-4), “ISG15- protein conjugation” (ISG15, UBA7 and UBE2E2, cor - rected p value 3.02 × 10-5) and “Homophilic cell adhesion via plasma membrane adhesion molecules” (AMIGO1, CDH15, CDH24, CDH6, FAT1, FAT2, PALLD, PCDHA9 and PLXNB3 were down-, while CDHR1 and NEC- TIN4 were upregulated, corrected p value 5.75 × 10-4). Upregulated genes were mapped in the specific network group “Cysteine metabolic process” (MPST , TST and VSIG2 were up- and SLC7A11 was downregulated, cor - rected p value 5.62 × 10-4). Nonspecific network groups characterized by equal proportions of mapped up- and downregulated genes were “Diseases associated with O-glycosylation of proteins” (ADAMTS17, ADAMTS5 and ADAMTSL2 were up-, while ADAMTSL1 , MUC13, MUC5B and THSD7A were downregulated, corrected p value 4.48 × 10-4) and “Retina homeostasis” (AZGP1, CDHR1 and NECTIN4 were up-, while ALPK3, ATP1B2, CDH15 and POTEJ were downregulated, corrected p value 6.42 × 10-4). The 33 identified enriched pathways are summarized in Additional file 5. Enriched pathways obtained by integration of identified DEGs and endometrial receptivity genes from the literature Only a set of 382 DEGs associated with the adenomyosis group that were identified by transcriptomics data com - parison of adenomyosis case and control samples dated to the receptive phase were used for integrative enrich - ment analyses with endometrial receptivity genes from the literature. Table 4 Top downregulated mRNAs and lncRNAs DEGs were insignificant after multiple testing correction of the p value. Abbreviation “FC” refers to fold change of expression levels ENTREZ ID HGNC symbol Long name Locus type logFC p value 169693 TMEM252 transmembrane protein 252 mRNA -2.1611 0.0331 5655 KLK10 kallikrein related peptidase 10 mRNA -1.7433 0.0010 5803 PTPRZ1 protein tyrosine phosphatase receptor type Z1 mRNA -1.7022 0.0120 727897 MUC5B mucin 5B, oligomeric mucus/gel-forming mRNA -1.6678 0.0421 79937 CNTNAP3 contactin associated protein like 3 mRNA -1.5039 0.0096 10752 CHL1 cell adhesion molecule L1 like mRNA -1.5027 0.0391 7103 TSPAN8 tetraspanin 8 mRNA -1.4535 0.0070 9723 SEMA3E semaphorin 3E mRNA -1.4216 0.0001 10964 IFI44 L interferon induced protein 44 like mRNA -1.3304 0.0416 5340 PLG plasminogen mRNA -1.3219 0.0499 145837 DRAIC downregulated RNA in cancer, inhibitor of cell invasion and migration lncRNA -1.9848 0.0392 100131825 CADM3-AS1 CADM3 antisense RNA 1 lncRNA -1.6197 0.0298 100506674 MRPS30-DT MRPS30 divergent transcript lncRNA -1.2230 0.0049 641364 SLC7A11-AS1 SLC7A11 antisense RNA 1 lncRNA -0.7862 0.0302 100506305 LINC00958 long intergenic non-protein coding RNA 958 lncRNA -0.6847 0.0448 100289410 MCF2 L-AS1 MCF2 L antisense RNA 1 lncRNA -0.6591 0.0414 386597 RNF144A-AS1 RNF144A antisense RNA 1 lncRNA -0.6095 0.0177 144481 SOCS2-AS1 SOCS2 antisense RNA 1 lncRNA -0.4608 0.0205 100134229 KDM7A-DT KDM7A divergent transcript lncRNA -0.4579 0.0344 Page 10 of 16Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 Integration of lists with 382, 151 and 173 genes associ - ated with adenomyosis, healthy uterus and endometrio - sis, respectively, provided 40 enriched pathways sorted in 11 network groups, which are presented in Fig.  3. According to the generated network, unique fingerprints of gynaecological pathologies on endometrial signa - tures were observed. The identified “Expression of IFN- induced genes” , “Negative regulation of viral process” and “Diseases associated with O-glycosylation of proteins” network groups were specific for the adenomyosis gene list, while “Interleukin-10 signalling” and “ ARC gene expression” were specific for the endometriosis gene list. In addition, nonspecific network groups, characterized by mapped genes originating from all 3 lists associated with gynaecological conditions, were identified, including “Extracellular matrix organization” , “Serine-type pepti - dase activity” , “Positive regulation of DNA-binding tran- scription factor activity” , “Cellular response to vascular endothelial growth factor stimulus” , “Response to cad - mium ion” and “Regulation of reproductive process” . This could indicate the interference of adenomyosis and endo- metriosis with molecular mechanisms required for nor - mal endometrial receptivity. The 40 identified enriched pathways are summarized in Additional file 6. Integration of lists with 424 (382 DEGs of the present sequencing experiment and 42 genes from the literature), 151 and 173 genes associated with adenomyosis, healthy uterus and endometriosis, respectively, provided 57 enriched pathways sorted in 18 network groups, which are presented in Fig.  4. Similar results were retrieved Fig. 2 Networks of enriched pathways and mapped genes associated with DEGs in the adenomyosis group. (a) Sorted 4 pathways in 2 network groups obtained by the enrichment analysis of 909 DEGs associated with adenomyosis group after comparing receptive, early- and late-receptive case and control samples; (b) Sorted 33 enriched pathways within 7 network groups identified from 382 DEGs associated with adenomyosis group after comparing only case and control samples dated to the receptive phase; (c) The enlargement of the connected network groups “Expression of IFN-induced genes” , “Response to interferon-alpha” and “ISG15-protein conjugation” presenting candidate pathways for future studies associated with altered endometrial receptivity in adenomyosis. Each set of DEGs was uploaded in the Cytoscape ClueGO app as two separate clusters, where upregulated genes were marked with violet and downregulated genes with green colour. Shape of nodes in networks attributed to ontology sources that were applied for enrichment analysis. Enriched pathways were sorted into network groups based on their common biological role Page 11 of 16 Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 with the integration of the adenomyosis gene list with 382 DEGs alone. However, some additional nonspecific network groups were identified, including “Interleukin-4 and Interleukin-13 signalling” , “Tumour necrosis fac - tor production” and “Sodium ion export across plasma membrane” . The 57 identified enriched pathways are summarized in Additional file 7.

Altered endometrial molecular mechanisms obstructing successful embryo implantation in women with adeno - myosis are poorly understood. The focus of the present transcriptomics analysis was to apply two novel molecu - lar approaches to identify gene expression differences in LH-timed endometrial samples between women with and without adenomyosis: genome-wide profiling using RNA-seq and accurate classification of endometrial receptivity as assessed by the molecular tool beREADY ®, measured from the same biopsy. Lists of DEGs associ - ated with the adenomyosis group that were identified by analysing RNA-seq datasets in the setting of the endome- trial dating results were applied for enrichment pathway analyses to predict their role in the context of endome - trial molecular organization. In addition, a set of 382 DEGs obtained after transcriptomics data comparison of confirmed receptive samples was used for further bioin - formatics analysis. They were integrated with the most extensive set of genes from the literature associated with endometrial receptivity in healthy uterus, endometrio - sis (model disease to study persistence of gynaecological pathology on endometrial molecular organization) and adenomyosis to predict mechanisms in which adenomyo- sis mediates an effect on endometrial receptivity. Recently, Devesa-Peiro et  al. [38] compared avail - able transcriptomics data and observed a greater effect of changing phases of the menstrual cycle on the Fig. 3 Integration of 382 adenomyosis-associated DEGs with endometrial receptivity genes associated with endometriosis and healthy uterus. In total, 40 pathways sorted into 11 network groups were obtained after enrichment analysis of integrated gene lists. The adenomyosis gene list (blue colour) included 382 DEGs associated with the adenomyosis group of the present RNA-seq analysis. The healthy uterus list (green colour) and endometriosis list (pink colour) contained 151 and 173 genes, respectively, which were associated with endometrial receptivity in the literature Page 12 of 16Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 endometrial transcriptome signature than the persis - tence of endometrial pathologies. Accurate dating of collected biopsies was highlighted as an important pro - cedure when identifying endometrial biomarkers asso - ciated with uterine abnormalities [38]. In addition, displacement of the temporal appearance of the WOI has been observed in some women [39, 40], meaning that the WOI could appear earlier or later in the luteal phase, as it is generally assumed that it is constant in all women [39, 41]. In view of these data, we utilized the novel molecular beREADY® tool [25], which reliably determines endome- trial dating on a transcriptomics platform, and machine- learning algorithms to assure homozygosity of LH-timed biopsies in the present study groups. Considering the

of endometrial receptivity testing, we excluded early- and late-receptive samples from the RNA-seq dataset to prevent the impact of early- and late-secretory phases associated with physiological advancement of endometrial maturation through the menstrual cycle, which could bias transcriptomics analysis associated with endometrial receptivity in adenomyosis. In that way, we identified 382 DEGs that we believe more accurately represent the effect of adenomyosis on the gene expres - sion signature of endometrial receptivity compared to 909 DEGs associated with the adenomyosis group, which were identified by comparing transcriptomics data of samples derived from receptive, early- and late-receptive phases. Enrichment analysis using 382 DEGs also provided a higher number of pathways tightly sorted in connected network groups compared to analysis of 909 DEGs, which were also more meaningful to relate with endo - metrial molecular biology (Fig.  2b). Namely, accord - ing to the results of the enrichment analysis of 382 Fig. 4 Integration of 424 adenomyosis-associated DEGs with endometrial receptivity genes associated with endometriosis and healthy uterus. In total, 57 enriched pathways sorted into 18 network groups were obtained after integration of adenomyosis list with 424 genes (382 DEGs of the present experiment and 42 genes from the literature), (blue colour), endometriosis list with 173 genes (pink colour) and healthy uterus list with 151 genes (green colour) Page 13 of 16 Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 DEGs, “Expression of IFN-induced genes” , “Response to interferon-alpha” and “ISG15-protein conjugation” were sorted as connected processes (Fig.  2c). Popovici et al. [42] associated increased expression levels of genes encoding chemotactic factors, inflammatory cytokines (including type I IFN-alpha/beta) and apoptosis-inducing agents with a role in the recruitment of lymphocytes and macrophages in human endometrial decidua [42]. IFNs, as reviewed by De Veer et al. [43], are a family of multi - functional cytokines that activate the expression of many genes with antiviral, antiproliferative or immunosuppres- sive effects. The signal transduction pathway of IFNs is initiated upon IFN binding to specific cell surface recep - tors. Downstream formed complexes of phosphorylated proteins and transcription factors bind to IFN-stimulated response elements (ISREs) at the promotor region of IFN-stimulated genes (ISGs) and initiate their transcrip - tion. There are more than 300 ISGs [43]. The ubiquitin- like protein ISG15 is a posttranscriptional modifier that can be in a process termed ISGylation covalently linked to hundreds of proteins. The role of ISG15 has been asso- ciated with cellular processes such as protein translation, cytoskeleton dynamics, exosome secretion, autophagy, genome stability and cancer; therefore, it presents a potential target for therapeutic strategies [44]. ISG15 can exert functions as an intracellular and secreted protein. Intracellular expression of ISG15, which is dependent on type I IFN-alpha/beta signalling, characterizes innate immune responses to viral and microbial pathogens. Its extracellular signalling can elicit secretion of cytokine type II IFN-gamma from lymphocytes [45]. Studies in mice suggested that ISG15 plays a role in the recruitment of uterine natural killer (uNK) cells during early gestation, where it is responsible for remodelling of spiral arteries to ensure a normal blood supply to the foetus and placenta throughout pregnancy [46]. The identified enriched path- ways related to the response to IFN signalling could indi - cate altered immune factors that have been associated with adenomyosis. Tremellen and Russell [47] associated an increased density of uNK cells and macrophages in the functional layer of late-secretory endometrium in women with severe adenomyosis experiencing implantation failures with a hostile immune environment that might interfere with successful embryo implantation [47]. In addition, Sotnikova et  al. [48] reported higher levels of secreted proinflammatory cytokines (IFN-gamma, IFN- alpha, tumour necrosis factor (TNF)-alpha and inter - leukin (IL)-1 beta) in supernatant samples of cultured mononuclear cells obtained from late-secretory endome - trium of women with adenomyosis when compared with healthy controls [48]. Another interesting enriched path - way from the 382 DEGs was related to cellular adhesion, whose importance in the process of embryo implantation has been described elsewhere [49]. The 382 Identified DEGs were also applied for the integration approach to repeat our previous enrichment pathway analysis [18] oriented to detect candidate path - ways of affected endometrial receptivity in adenomyosis. Integrative enrichment analysis using the adenomyosis gene list with 382 DEGs only provided candidate path - ways associated with endometrial receptivity establish - ment (e.g., “Extracellular matrix organization” “Cellular response to vascular endothelial growth factor stimulus” and “Regulation of reproductive process”) that could be dysregulated in adenomyosis as well as in endome - triosis, which is in agreement with the literature [17, 50–54]. The identified specific network group “Expres - sion of IFN-induced genes” persisted as a unique effect of adenomyosis on endometrial molecular background after enrichment analysis using integrated gene lists. Enriched pathways related to activity-regulated cytoskeletal (ARC) gene expression were specific to the endometriosis gene list, which was used as a model to study the effect of endometrial-associated disorders. ARC is an immediate early gene involved in signal transduction. Its transcrip - tion is induced by various signalling cascades, includ - ing mitogen-activated protein kinases (MAPKs) and extracellular signal-regulated kinases (ERKs) [55], which have already been associated with endometrial receptiv - ity defects in endometriosis [56]. Integrative enrichment analysis using the adenomyosis gene list with 424 genes provided additional candidate pathways to be associ - ated with altered cytokine responses in adenomyosis and endometriosis, including “interleukin-4 and interleu - kin-13 signalling” , which was also identified in our previ- ous study [18] and could be attributable to the Reactome pathway database being used as an ontology source in both studies, “regulation of TNF superfamily cytokine production” and “interleukin-10 signalling” . Altered expression levels of some cytokines in the endometrium during WOI have been observed in women with adeno - myosis after COS [57] and in women with endometriosis [58, 59]. Decidualization of endometrial stromal cells is characterized by a changing endometrial inflammatory environment shown as a transition from a proinflamma - tory to an anti-inflammatory response [60, 61]. This tran- sition has been associated with balancing endometrial receptivity versus selectively accepting only high-quality embryos [61]. Dysregulated balance has been associated with the implantation of poor-quality embryos leading to miscarriage [62]. It could be that enriched pathways asso- ciated with the expression of IFN-induced genes indicate dysregulated endometrial selectively, which may explain the observed higher incidence of early pregnancy loss in women with adenomyosis [6, 20]. However, further Page 14 of 16Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 studies are needed to verify this hypothesis. The identi - fication of robust pathways could lead to the extension of current gene sets for endometrial receptivity examination presented in a growing number of commercial molecular tests [25, 63–65] that would be specific for women with adenomyosis. Accurate endometrial receptivity exami - nation in this group of infertile patients could better verify whether endometrial-associated factor is a source of recurrent implantation failures that prolong infertil - ity treatments [66]. Furthermore, endometrial changes in women with adenomyosis could provide not only the relationship between pathophysiological mechanisms of adenomyosis development [17, 67, 68] but also the patho- genesis of the malignant transformation [69]. Recently, it was reported that endometrial carcinoma could co-exist or arise from adenomyosis which may be important fac - tor in survival outcomes of the patient [70, 71]. A limitation of our study is the relatively small sam - ple size, which prevents definitive conclusions regarding the impact of adenomyosis on the endometrial tran - scriptome [72]. The results could also differ if the con - trol group is composed of women with proven fertility. Another limitation of the present study is that the diag - nosis of adenomyosis could only be made noninvasively by imaging, since definitive histopathological diagno - sis can only be made after hysterectomy. In genome- wide studies focusing on pathophysiological aspects of adenomyosis, the diagnosis can be based on histological examination of specimens after hysterectomy [17, 67, 68]. However, this is only possible retrospectively and is irrelevant in women who wish to preserve their fer - tility. In fertility-oriented transcriptomics studies [20] or studies including endometriosis [73], a diagnosis of adenomyosis was noninvasive. The diagnosis of adeno - myosis by ultrasound is challenging, and there are no uniform ultrasonographic criteria for the diagnosis [74]. In the present study, TVUS of the uterus and pelvic cav - ity was performed by an experienced sonographer prior to each endometrial biopsy to confirm sonographic evidence of adenomyosis and to exclude other pelvic pathologies.

In this study, we focused on the molecular background of infertility-related adenomyosis based on our research and the available literature. We applied accurate endo - metrial receptivity classification of retrieved endometrial samples LH-timed to the expected WOI to avoid men - strual cycle bias in downstream transcriptomics analy - sis. The 382 DEGs identified in the adenomyosis group using the RNA-seq dataset of only confirmed receptive endometrial samples resulted in 33 enriched pathways further projected in the network from which “Expres - sion of IFN-induced genes” , “Response to interferon- alpha” and “ISG15-protein conjugation” were highlighted as connected processes. Additional integration of 382 DEGs with candidate genes associated with endome - trial receptivity in healthy uterus, endometriosis and adenomyosis based on a literature review revealed that cytokine signalling impairments in endometrial patholo - gies could interfere with mechanisms of endometrial receptivity. According to our results, an altered response to IFN signalling is suggested as a candidate mechanism of impaired uterine receptivity in adenomyosis that needs to be further studied in a larger sample size. Abbreviations ART : Assisted reproductive technique; DEG: Differentially expressed gene; GO: Gene ontology; IFN: Interferon; LH: Luteinizing hormone; RNA-seq: RNA sequencing; TVUS: Transvaginal ultrasound ; WOI: Window of implantation. Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s12958- 021- 00871-5. Additional file 1. Additional file 2. Additional file 3. Additional file 4. Additional file 5. Additional file 6. Additional file 7.

The authors are grateful to the patients for their participation in this study. We thank Jasna Muršič for her assistance in the recruitment of patients. Authors’ contributions EP , TK, BK and JK designed research. EP , MG, JK, TK analysed the data. EP and JK wrote the manuscript. TK, JK, BK, and UP revised the manuscript. All authors have read and approved the final manuscript. Funding This work was a part of research programmes P3–0327 and P4–0220, funded by the Slovenian Research Agency (ARRS). Availability of data and materials All data generated and analysed during this study are included in this article [and its supplementary information file]. The RNA-seq data presented in this study are deposited in the GEO database with accession number GSE185392. Declarations Ethics approval and consent to participate The study was approved by the Slovenian National Medical Ethics Committee (0120–259/2018/16). Each patient signed an informed consent form before being involved in the study. Page 15 of 16 Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 Competing interests The authors declare that they have no competing interests. Author details 1 Department of Reproductive Medicine and Gynaecological Endocrinology, University Medical Centre Maribor, 2000 Maribor, Slovenia. 2 Department of Animal Science, Biotechnical Faculty, University of Ljubljana, 1230 Domžale, Slovenia. 3 Centre for Human Molecular Genetics and Pharmacogenomics, Faculty of Medicine, University of Maribor, 2000 Maribor, Slovenia. 4 Labora- tory for Biochemistry, Molecular Biology and Genomics, Faculty of Chemistry and Chemical Engineering, University of Maribor, 2000 Maribor, Slovenia. 5 Department of Gynaecology, University Medical Centre Maribor, 2000 Mari- bor, Slovenia. Received: 23 August 2021 Accepted: 6 October 2021

1. Andres MP , Borrelli GM, Ribeiro J, Baracat EC, Abrão MS, Kho RM. Trans- vaginal ultrasound for the diagnosis of adenomyosis: systematic review and meta-analysis. J Minim Invasive Gynecol [Internet]. Elsevier Inc.; 2018;25:257–64. Available from: https:// doi. org/ 10. 1016/j. jmig. 2017. 08. 653 2. Pinzauti S, Lazzeri L, Tosti C, Centini G, Orlandini C, Luisi S, et al. Trans- vaginal sonographic features of diffuse adenomyosis in 18-30-year-old nulligravid women without endometriosis: association with symptoms. Ultrasound Obstet Gynecol. 2015;46:730–6. 3. Puente JM, Fabris A, Patel J, Patel A, Cerrillo M, Requena A, et al. Adeno- myosis in infertile women: prevalence and the role of 3D ultrasound as a marker of severity of the disease. Reprod Biol Endocrinol [Internet]. Reproductive Biology and Endocrinology; 2016;14:60. Available from: https:// doi. org/ 10. 1186/ s12958- 016- 0185-6 4. Sharma S, Bathwal S, Agarwal N, Chattopadhyay R, Saha I, Chakravarty B. Does presence of adenomyosis affect reproductive outcome in IVF cycles? A retrospective analysis of 973 patients. Reprod Biomed Online. Elsevier Ltd; 2019;38:13–21. 5. Thalluri V, Tremellen KP . Ultrasound diagnosed adenomyosis has a nega- tive impact on successful implantation following GnRH antagonist IVF treatment. Hum Reprod. 2012;27:3487–92. 6. Stanekova V, Woodman RJ, Tremellen K. The rate of euploid miscar- riage is increased in the setting of adenomyosis. Hum Reprod Open. 2019;2019:1–8. 7. Mehasseb MK, Bell SC, Pringle JH, Habiba MA. Uterine adenomyosis is associated with ultrastructural features of altered contractility in the inner myometrium. Fertil Steril [Internet]. Elsevier Ltd; 2010;93:2130–6. Avail- able from: https:// doi. org/ 10. 1016/j. fertn stert. 2009. 01. 097 8. Harada T, Khine Y, Kaponis A, Nikellis T, Decavalas G, F T. The impact od adenomyosis on women’s fertility. Obstet Gynecol Surv. 2016;71:557–68. 9. Ota H, Igarashi S, Hatazawa J, Tanaka T. Endothelial nitric oxide synthase in the endometrium during the menstrual cycle in patients with endome- triosis and adenomyosis. Fertil Steril. 1998;69:303–8. 10. Ota H, Igarashi S, Hatazaw J, Tanaka T. Immunohistochemical assessment of superoxide dismutase expression in the endometrium in endometrio- sis and adenomyosis. Fertil Steril. 1999;72:129–34. 11. Ota H, Igarashi S, Sasaki M, Tanaka T. Xanthine oxidase in eutopic and ectopic endometrium in endometriosis and adenomyosis. Hum Reprod. 2001;75:785–90. 12. Ota H, Tanaka T. Stromal vascularization in the endometrium during adenomyosis. Microsc Res Tech. 2003;449:445–9. 13. Goteri G, Lucarini G, Montik N, Zizzi A, Stramazzotti D, Fabris G, et al. Expression of vascular endothelial growth factor (VEGF), hypoxia induc- ible factor-1α (HIF-1α), and microvessel density in endometrial tissue in women with adenomyosis. Int J Gynecol Pathol. 2009;28:157–63. 14. Peng Y, Jin Z, Liu H, Xu C. Impaired decidualization of human endometrial stromal cells from women with adenomyosis. Biol Reprod. 2021;ioab017. 15. Wang J, Deng X, Yang Y, Yang X, Kong B, Chao L. Expression of GRIM- 19 in adenomyosis and its possible role in pathogenesis. Fertil Steril. 2016;105:1093–101. 16. Campo S, Campo V, Benagiano G. Adenomyosis and infertility. Reprod Biomed Online. Reproductive Healthcare Ltd.; 2012;24:35–46. 17. Herndon CN, Aghajanova L, Bayalan S, Erikson D, Barragan F, Goldfien G, et al. Global transcriptome abnormalities of the eutopic endometrium from women with adenomyosis. Reprod Sci. 2016;23:1289–303. 18. Prašnikar E, Kunej T, Repnik K, Potočnik U, Knez J, Kovačič B. Determining the molecular background of endometrial receptivity in adenomyosis. Biomolecules. 2020;10:1–25. 19. Prašnikar E, Knez J, Kovačič B, Kunej T. Molecular signature of eutopic endometrium in endometriosis based on the multi-omics integrative synthesis. J Assist Reprod Genet. 2020;37:1593–611. 20. Martinez-Conejero JA, Morgan M, Montesinos M, Fortu S, Meseguer M, Ph D, et al. Adenomyosis does not affect implantation, but is associated with miscarriage in patients undergoing oocyte donation. Fertil Steril. 2011;96:943–51. 21. Zhao S, Fung-Leung WP , Bittner A, Ngo K, Liu X. Comparison of RNA-Seq and microarray in transcriptome profiling of activated T cells. PLoS One. 2014;9:e78644. 22. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for tran- scriptomics. Nat Rev Genet. 2009;10:57–63. 23. Naftalin J, Hoo W, Pateman K, Mavrelos D, Holland T, Jurkovic D. How common is adenomyosis? A prospective study of prevalence using trans- vaginal ultrasound in a gynaecology clinic. Hum Reprod. 2012;27:3432–9. 24. Hulka CA, Hall DA, McCarthy K, Simeone J. Sonographic findings in patients with adenomyosis: can sonography assist in predicting extent of disease? Am J Roentgenol. 2002;179:379–83. 25. Saare M, Laisk T, Teder H, Paluoja P , Palta P , Koel M, et al. A molecular tool for menstrual cycle phase dating of endometrial samples in endometrio- sis transcriptome studies. Biol Reprod. 2019;101:1–3. 26. Teder H, Koel M, Paluoja P , Jatsenko T, Rekker K, Laisk-Podar T, et al. TAC- seq: targeted DNA and RNA sequencing for precise biomarker molecule counting. NPJ Genomic Med. 2018;3. 27. Altmäe S, Koel M, Võsa U, Adler P , Suhorutšenko M, Laisk-Podar T, et al. Meta-signature of human endometrial receptivity: a meta-analysis and validation study of transcriptomic biomarkers. Sci Rep [Inter- net]. 2017;7:10077. Available from: http:// www. nature. com/ artic les/ s41598- 017- 10098-3 28. Liao Y, Smyth GK, Shi W. The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res. Oxford University Press; 2019;47. 29. Liao Y, Smyth GK, Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 2013;41:e108. 30. Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformat- ics. 2014;30:923–30. 31. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioin- formatics. 2010;26:139–40. 32. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25. 33. Law CW, Chen Y, Shi W, Smyth GK. Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15:R29. 34. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray stud- ies. Nucleic Acids Res. 2015;43:e47. 35. Bindea G, Mlecnik B, Hackl H, Charoentong P , Tosolini M, Kirilovsky A, et al. ClueGO: A Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25:1091–3. 36. Bindea G, Galon J, Mlecnik B. CluePedia Cytoscape plugin: Pathway insights using integrated experimental and in silico data. Bioinformatics. 2013;29:661–3. 37. Shannon P , Markiel A, Ozier O, Baliga SN, Wang TJ, Ramage D, et al. Cytoscape: a software environment for integrated models of biomo- lecular interaction networks. Genome Res [Internet]. 2003;13:2498–504. Available from: http:// ci. nii. ac. jp/ naid/ 11000 19104 81/ 38. Devesa-Peiro A, Sebastian-Leon P , Pellicer A, Diaz-Gimeno P . Guidelines for biomarker discovery in endometrium: correcting for menstrual cycle bias reveals new genes associated with uterine disorders. Mol Hum Reprod. 2021;27:1–19. Page 16 of 16Prašnikar et al. Reproductive Biology and Endocrinology (2022) 20:2 39. Ruiz-Alonso M, Blese D, Diaz-Gimeno P , Gomez E, Fernandez-Sanchez M, Carranza F, et al. The endometrial receptivity array for diagnosis and personalized embryo transfer as a treatment for patients with repeated implantation failure. Fertil Steril. 2013;100:818–24. 40. Mahajan N, Kaur S, Alonso MR. Window of implantation is significantly displaced in patients with adenomyosis with previous implantation failure as determined by endometrial receptivity assay. J Hum Reprod Sci. 2018;11:353–8. 41. Tan J, Kan A, Hitkari J, Taylor B, Tallon N, Warraich G, et al. The role of the endometrial receptivity array (ERA) in patients who have failed euploid embryo transfers. J Assist Reprod Genet. Journal of Assisted Reproduction and Genetics; 2018;35:683–92. 42. Popovici RM, Kao LC, Giudice CG. Discovery of new inducible genes in in vitro decidualized human endometrial stromal cells using microarray technology. Endocrinology. 2000;141:3510–3. 43. de Veer MJ, Holko M, Frevel M, Walker E, Der S, Paranjape JM, et al. Functional classification of interferon-stimulated genes identified using microarrays. J Leukoc Biol [Internet]. 2001;69:912–20. Available from: http:// www. ncbi. nlm. nih. gov/ pubmed/ 11404 376 44. Jiménez Fernández D, Hess S, Knobeloch KP . Strategies to target ISG15 and USP18 toward therapeutic applications. Front Chem. 2020;21:923. 45. Young Hee Choi and A-MY. Extracellular ISG15 signals cytokine secretion through the LFA-1 integrin receptor. Mollecular cell. 2017;68:581-590.e5. 46. Austin KJ, Bany BM, Belden EL, Rempel LA, Cross JC, Hansen TR. Interferon-stimulated gene-15 (Isg15) expression is up-regulated in the mouse uterus in response to the implanting conceptus. Endocrinology. 2003;144:3107–13. 47. Tremellen KP , Russell P . The distribution of immune cells and mac- rophages in the endometrium of women with recurrent reproductive failure. II: adenomyosis and macrophages. J Reprod Immunol [Internet]. Elsevier Ireland Ltd; 2012;93:58–63. Available from: https:// doi. org/ 10. 1016/j. jri. 2011. 12. 001 48. Sotnikova N, Antsiferova I, Malyshkina A. Cytokine network of eutopic and ectopic endometrium in women with adenomyosis. Am J Reprod Immunol. 2002;47:251–5. 49. Achache H, Revel A. Endometrial receptivity markers, the journey to suc- cessful embryo implantation. Hum Reprod Update. 2006;12:731–46. 50. Matsuzaki S, Maleysson E, Darcha C. Analysis of matrix metallopro- teinase-7 expression in eutopic and ectopic endometrium samples from patients with different forms of endometriosis. Hum Reprod. 2010;25:742–50. 51. Pino M, Galleguillos C, Torres M, Sovino H, Fuentes A, Boric MA, et al. Association between MMP1 and MMP9 activities and ICAM1 cleav- age induced by tumor necrosis factor in stromal cell cultures from eutopic endometria of women with endometriosis. Reproduction. 2009;138:837–47. 52. Okada H, Tsuzuki T, Shindoh H, Nishigaki A, Yasuda K, Kanzaki H. Regula- tion of decidualization and angiogenesis in the human endometrium: mini review. J Obstet Gynaecol Res. 2014;40:1180–7. 53. Li T, Li Y-G, Pu D-M. Matrix metalloproteinase-2 and -9 expression cor- related with angiogenesis in human adenomyosis. Gynecol Obstet Invest. 2006;62:229–35. 54. Di Carlo C, Bonifacio M, Tommaselli GA, Bifulco G, Guerra G, Nappi C. Met- alloproteinases, vascular endothelial growth factor, and angiopoietin 1 and 2 in eutopic and ectopic endometrium. Fertil Steril [Internet]. Elsevier Ltd; 2009;91:2315–23. Available from: https:// doi. org/ 10. 1016/j. fertn stert. 2008. 03. 079 55. Epstein I, Finkbeiner S. The Arc of cognition: signaling cascades regulating Arc and implications for cognitive function and disease. Semin Cell Dev Biol. 2018;77:63–72. 56. Lessey BA, Kim JJ. Endometrial receptivity in eutopic endometrium of women with endometriosis it is affected, let me show you why. Fertil Steril. 2017;108:19–27. 57. Zhihong N, Yun F, Pinggui Z, Sulian Z, Zhang A. Cytokine profiling in the eutopic endometrium of adenomyosis during the implantation window after ovarian stimulation. Reprod Sci. 2016;23:124–33. 58. Kharfi A, Labelle Y, Mailloux J, Akoum A. Deficient expression of tumor necrosis factor receptor type 2 in the endometrium of women with endometriosis. Am J Reprod Immunol. 2003;50:33–40. 59. Ulukus M, Ulukus EC, Seval Y, Zheng W, Arici A. Expression of interleukin-8 receptors in endometriosis. Hum Reprod [Internet]. 2005;20:794–801. Available from: https:// acade mic. oup. com/ humrep/ artic le- lookup/ doi/https:// doi. org/ 10. 1093/ humrep/ deh675 60. Wilczyński JR. Th1/Th2 cytokines balance - yin and yang of reproductive immunology. Eur J Obstet Gynecol Reprod Biol. 2005;122:136–43. 61. Macklon NS, Brosens JJ. The human endometrium as a sensor of embryo quality. Biol Reprod. 2014;91:1–8. 62. Quenby S, Vince G, Farquharson R, Aplin J. Recurrent miscarriage: a defect in nature’s quality control? Hum Reprod. 2002;17:1959–63. 63. Díaz-Gimeno P , Horcajadas JA, Martinez-Conejero JA, Esteban FJ, Alama P , Pellicer A, et al. A genomic diagnostic tool for human endometrial recep- tivity based on the transcriptomic signature. Fertil Steril. 2011;95:50–60. 64. Enciso M, Carrascosa JP , Sarasa J. Development of a new comprehensive and reliable endometrial receptivity map (ER Map / ER Grade) based on RT- qPCR gene expression analysis. Hum Reprod. 2018;33:220–8. 65. Haouzi D, Entezami F, Torre A, Innocenti C, Antoine Y, Mauries C, et al. Customized frozen embryo transfer after identification of the receptivity window with a transcriptomic approach improves the implantation and live birth rates in patients with repeated implantation failure. Reprod Sci. Reproductive Sciences; 2021;28:69–78. 66. Tremellen K, Russell P . Adenomyosis is a potential cause of recurrent implantation failure during IVF treatment. Aust New Zeal J Obstet Gynae- col. 2011;51:280–3. 67. Jiang JF, Sun AJ, Xue W, Deng Y, Wang YF. Aberrantly expressed long noncoding RNAs in the eutopic endometria of patients with uterine adenomyosis. Eur J Obstet Gynecol Reprod Biol [Internet]. Elsevier Ireland Ltd; 2016;199:32–7. Available from: https:// doi. org/ 10. 1016/j. ejogrb. 2016. 01. 033 68. Xiang Y, Sun Y, Yang B, Yang Y, Zhang Y, Yu T, et al. Transcriptome sequenc- ing of adenomyosis eutopic endometrium: a new insight into its patho- physiology. J Cell Mol Med. 2019;23:8381–91. 69. Laganà AS, Scioscia M. Endometrial cancer in women with adenomyosis: an underestimated risk? Int J Fertil Steril. 2020;14:260–1. 70. Chao X, Wu M, Ma S, Tan X, Zhong S, Bi Y, et al. The clinicopathological characteristics and survival outcomes of endometrial carcinoma coexist- ing with or arising in adenomyosis: a pilot study. Sci Rep. 2020;10:1–9. 71. Scioscia M, Noventa M, Laganà AS. Abnormal uterine bleeding and the risk of endometrial cancer: can subendometrial vascular ultrasound be of help to discriminate cancer from adenomyosis? Am J Obstet Gynecol [Internet]. Elsevier Inc.; 2020;223:605–6. Available from: https:// doi. org/ 10. 1016/j. ajog. 2020. 05. 049 72. Vercellini P , Parazzini F, Oldani S, Panazza S, Bramante T, Crosignani PG. Adenomyosis at hysterectomy: a study on frequency distribution and patient characteristics. Hum Reprod. 1995;10:1160–2. 73. Dior UP , Nisbet D, Fung JN, Foster G, Healey M, Montgomery GW, et al. The association of sonographic evidence of adenomyosis with severe endometriosis and with gene expression in eutopic endometrium. J Minim Invasive Gynecol [Internet]. Elsevier Inc.; 2018;26:941–8. Available from: https://linkinghub.elsevier.com/ retrieve/pii/S1553465018312512%0Ahttp://www.ncbi.nlm.nih.gov/ pubmed/30273686 74. Chapron C, Vannuccini S, Santulli P , Abrão MS, Carmona F, Fraser IS, et al. Diagnosing adenomyosis: an integrated clinical and imaging approach. Hum Reprod Update. 2020;26:392–411. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub- lished maps and institutional affiliations.