The study was approved by the Research Ethics Committee of the University of Tartu, Estonia (approval No 333/T-6), and all the study participants have provided written informed consent. The tissue samples were collected from patients undergoing surgery at the Tartu University Hospital (Tartu, Estonia) and with confirmed endometriosis in accordance with the revised American Society for Reproductive Medicine classification system (26). Women undergoing laparoscopy due to infertility or pelvic pain, and confirmed endometriosis-free, were enrolled as controls (women without endometriosis). Peritoneal lesions were removed by laparoscopic surgery, and the collection of endometrial biopsy samples from women with and without endometriosis was done with an endometrial suction catheter (Pipelle, Laboratoire CCD). The enrolled women had not received any hormonal therapy for at least 3 months before the sample collection. Table 1 provides detailed information about the patients enrolled in this study. The collected tissues were cut into two portions followed by endometriotic tissue histological evaluation and confirmation, or cryopreservation in Dulbecco’s Modified Eagle’s Medium (DMEM, Gibco, Thermo Fisher Scientific), with 30% fetal bovine serum (FBS, Biowest) and 7.5% dimethyl sulfoxide Hybri-Max (Sigma-Aldrich). The cryopreserved samples were stored in liquid nitrogen until further use.

Primary endometrial SCs were isolated from EuE, and EcE (peritoneal lesions) as previously described (3). EuESCs and EcESCs were cultured up to three passages in phenol red-free DMEM (Gibco) supplemented with 10% charcoal-stripped FBS (Gibco) with penicillin/streptomycin/amphotericin B (100 U/mL, 100 μg/mL, and 0.25 μg/mL, respectively) at 37°C under 20% O₂ and 5% CO₂. The cells were cryopreserved (1 vial from 1 confluent 10 cm Petri dish) and stored in liquid nitrogen until further use for functional assays.

The cells were thawed, seeded onto 6-well plates (working volume of 2 mL per well, approximately 170,000 cells per well) and cultured in phenol red-free DMEM with 10% charcoal-stripped FBS and penicillin/streptomycin/amphotericin B at 37°C under 20% O₂ and 5% CO₂ for 24 hours. Next, the medium was exchanged (working volume of 3 mL per well) and half of the replicates were cultured at 20% oxygen (normoxia) and the other half at 1% oxygen (hypoxia) for 48 hours. The hypoxic gas mixture (1% O₂, 5% CO₂, 94% N₂) was supplied by the gas controller into the incubator of Cytation 5 multi-mode reader system equilibrated at 37°C (BioTek).

After incubation, the medium was removed (1 mL aliquots of spent medium were frozen and subsequently used for the enzyme-linked immunosorbent assay, ELISA) and the wells were rinsed with PBS supplemented with Ca²⁺ and Mg²⁺ (Corning). The cells were then lysed either at room temperature using RLT buffer (Qiagen; for RNA extraction) or on ice using the Triton X-containing HEPES buffer (Calbiochem, Ferak; for protein kinase assay).

The fresh lysates prepared as previously described were used for protein kinase assay (20). For each lysate, total protein concentration was measured using Bradford assay (Thermo Scientific) according to the manufacturer’s instructions. Bradford assay was carried out in 96-well clear flat bottom plates (Nunc™ 269620) and absorbance at 590 nm was measured using a PHERAstar multi-mode reader (BMG Labtech). Within each experiment, the total protein concentration of all collected lysates was equalized by suitable dilution with kinase assay buffer as previously described (20).

The protein kinase assay was performed according to the previously published protocol using in-house synthesized photoluminescent probes that gain long-lifetime luminescence in complex with basophilic PKs (probe ARC-1139) or acidophilic PKs (probe ARC-1530) and the following displacing compounds: H89 (Biaffin) in case of PKAc; Y-27632 (Cayman Chemical) in case of ROCK; GSK690693 (Tocris) in case of Akt; and CX-4945 (Synkinase) in case of CK2 (20, 27, 28). The assay was carried out in 384-well round-bottom low volume plates (Corning 4514) and the time-delayed photoluminescence measurements were performed using PHERAstar multi-mode reader (BMG Labtech) with the following parameters: excitation 337 (300–360) nm, emission 675 (50) nm for ARC-1139 probe or 590 (50) nm for ARC-1530 probe, 100 flashes per well, integration start 60 μs, integration time 400 μs.

The total RNA was isolated from the SC culture lysates using RNeasy mini kit (Qiagen) following the manufacturer’s instructions. RNA samples were treated with DNase to remove DNA content with Ambion^® DNA-free™ DNase Treatment and Removal Reagents (ThermoFisher Scientific). RNA concentrations were measured with a Qubit Flex Fluorometer (Invitrogen) following the protocol from the Qubit High Sensitivity RNA assay kit.

The DNA libraries were constructed following the Smart-seq2 protocol described previously (29) with 10 ng of each RNA sample, then enzymatically fragmented and tagged using Nextera XT kit (Illumina Inc) and IDT^® for Illumina^® DNA Unique Dual Index barcodes. The purification of the final amplified cDNA libraries was performed using AMPure XP beads (1:1 ratio for the sample vs beads). The quantification of cDNA libraries was performed using the Qubit 1X HS DNA assay kit (Invitrogen), followed by the assessment of library quality using a high-sensitivity DNA chip (Agilent) on the 2100 Bioanalyzer system (Agilent). The pooled libraries were subjected to next-generation sequencing on the Illumina NovaSeq 6000 sequencing platform using a 2 x 150 bp read length (Novogene, UK).

The raw RNA-seq reads were subject to quality control assessment using the FastQC program (version 0.11.8) (30). Next, a comprehensive report based on the FASTQC results was generated using the MultiQC program (31) to evaluate the sequencing data quality. The fastp program (32) was then utilized with default parameters to remove adaptor sequences and perform read trimming. The resulting trimmed reads were aligned to the human reference genome GRCh38 using STAR 2.7.5a (33), accounting for paired-end reads. Gene counts were derived from the aligned reads utilizing the featureCounts program (34) with default parameters, also considering paired-end reads. On average, mRNA-seq generated 5.2 million reads per sample, ranging from 3.5 million to 7.2 million reads.

Principal component analysis (PCA) was performed using the plotPCA function available in the DESeq2 R package (35). Prior to performing PCA, the raw count data underwent variance stabilizing transformation (VST) to alleviate potential biases and improve clustering accuracy. For differential gene expression analysis, the DESeq2 R package (35) was utilized. We applied gene filtering to include only those genes with at least ten counts across all samples in at least one of the experimental groups. Next, pairwise differential gene expression analysis was conducted between experimental groups indicated in Table 2 . The resulting p-values underwent adjustment using Benjamini and Hochberg’s (BH) method to control the false discovery rate (FDR). Genes exhibiting an absolute log₂ fold change |log₂FC| > 0.5 and an FDR-adjusted P-value (P_adj) < 0.05 were considered differentially expressed.

Pathway analysis was performed using the g:Profiler web tool (36). Specifically, the g:GOSt function within g:Profiler facilitated gene-set enrichment analysis, with significance determined using the tailored significance threshold g:SCS. KEGG and GO: BP plots were made using the online tool ShinyGO 0.80 (37–39). Venn diagrams were made using the jvenn online tool (40). Only the genes with nomenclature names were included in the list used for creating the Venn diagram.

The amount of secreted TGFBI protein was quantified in spent cell culture media using the EHTGFBI (BIGH3) human ELISA kit (Invitrogen) according to the manufacturer’s instructions. Absorbance at 450 nm was measured using a Cytation 5 multi-mode reader (BioTek).

cDNAs were synthesized by RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific). The qRT-PCR analysis for transforming growth factor beta induced (TGFBI) expression was performed using 5x HOT FIREPol EvaGreen qPCR Mix Plus (ROX) (Solis BioDyne) with the primers (F: 5´-TAACGGCCAGTACACGCTT-3´; R: 5´-TGTTCAGCAGGTCTCTCAGG-3´). Gene expression data was normalized to the SDHA housekeeping gene (F: 5’-TGGGAACAAGAGGGCATCTG-3’; R: 5’-CCACCACTGCATCAAATTCATG-3’) as described previously (41).

To localize TGFBI protein in EuE and EcE, IHC analysis of tissue samples was performed. The tissue specimens were fixed in 10% neutral buffered formalin, dehydrated and embedded in paraffin at TUH’s Pathology Service. The prepared 2.5 µm tissue sections were deparaffinized in xylene followed by rehydration with a series of washes (2x, each for 5 min) in ethanol (96%, 90%, 80%, 70%, 60%, 50%) and MQ water. Next, tissue sections were incubated for 20 min at 98°C in 10 mM sodium-citrate buffer (pH = 6.0) for antigen retrieval and 5% normal goat serum in 1% BSA/PBS was used for 60 min to block nonspecific binding. The following steps were performed according to the Master Polymer Plus Detection System (Peroxidase) (Incl. 3,3’-Diaminobenzidine (DAB) Chromogen) Kit (Master Diagnostica, Spain) guidelines. All the washes were performed with 1x TBS (Fisher Bioreagents, USA), 3 times, each for 5 min. To identify TGFBI protein, tissue sections were incubated with rabbit anti-TGFBI antibody (Invitrogen, JM24-53, catalogue N MA5-32736, dilution 1:100) overnight at 4°C, and as negative controls, the tissue sections were incubated with 1% BSA/PBS. For nuclear counterstain, Mayer hematoxylin solution 1:4 was used. Lastly, the tissue sections were dehydrated in the opposite order as for the rehydration procedure, by incubating in ethanol and xylene. A mounting medium (Leica CV mount, Leica Biosystems, USA) was applied to the tissue sections and covered by coverslips. The tissue sections were then scanned with a Leica SCN400 slide scanner (Leica Biosystems, USA) using a 20x objective.

Semi-quantitative analysis of scanned sections was performed using the ImageJ Fiji package (v. 1.54k) (42). Briefly, DAB signal intensity was measured in the endometrial stroma, excluding luminal and glandular epithelium. The relative DAB intensity was calculated using the formula: f = 255 - i, where “f” is relative DAB intensity and “i” is mean DAB intensity obtained from the software ranging from 0 (deep brown, highest expression) to 255 (total white) (43).

PK assay: The initial normalization of signals was carried out according to the kinase activity profile observed in EuESCs under normoxic conditions (set to 100% for each PK of interest). Due to a large spread of data points observed for EcESCs, the second round of normalization was then carried out where EcESCs were normalized separately from EuESCs (cells treated under normoxia set to 100% for each PK of interest). The pairwise comparisons between the PK activity profiles in different oxygenation conditions or eutopic vs ectopic cell lysates (N = 5) were carried out using the unpaired two-tailed t-test with Welch’s correction and the significance of comparisons is indicated as follows: ** corresponds to P ≤ 0.01, * corresponds to P ≤ 0.05.

ELISA assay: For each sample of spent medium, the calculated TGFBI concentration was normalized to the total protein content measured for the corresponding cell lysate (i.e., homogenate of cells growing in the well from which the medium sample was initially collected) as described in Section 2.4. An additional normalization was carried out where for samples obtained from the same patient, the relative TGFBI amount was normalized according to the signal measured for the EuESCs under normoxia (set to 100%) or for the EuESCs or EcESCs incubated in normoxia (= 100%) separately. The pairwise comparisons between the oxygenation conditions or eutopic vs ectopic cells (N = 5) were carried out using the unpaired two-tailed t-test with Welch’s correction and the significance of comparisons is indicated as follows: *** corresponds to P ≤ 0.001, ** corresponds to P ≤ 0.01.

qRT-PCR: The comparison between TGFBI mRNA expression levels (ΔCt values) in proliferative vs secretory phase endometrium from women with (N = 45) and without (N = 24) endometriosis was performed using the unpaired two-tailed t-test with Welch’s correction and the significance of comparisons is indicated as follows: **** corresponds to P < 0.0001.

IHC: Statistical significance was assessed using the Wilcoxon Mann-Whitney test, with *corresponding to P ≤ 0.05.

For statistical analysis of data from ELISA, PK assay, qRT-PCR, and IHC, GraphPad Prism 6 (San Diego) was used as reported previously (20).

We started our study by examining the activity profiles of PKs that have been previously shown to contribute to the decidualization of normal endometrium – with PKA, ROCK and Akt/PKB featuring increased activity and CK2 featuring decreased activity in the ESCs decidualized in vitro or in situ. Our current goal was to establish whether the activity of these kinases can be altered depending on the origin of cells (eutopic vs ectopic) or short-term changes in oxygenation conditions (20% vs 1% oxygen). The activities of four PK families – PKAc, ROCK, Akt, and CK2 – were measured in fresh lysates prepared directly after incubation of EuESCs and EcESCs under different oxygenation conditions. To compensate for the differences in number of cells following the 48-h incubation, the lysates obtained from the same patients’ samples were diluted to the same total protein content. For each individual patient, the relative PK activity was normalized based on the PK activity established for EuESCs in normoxia, which was set to 100%.

The pooled results (N = 5) are summarized in Figure 1A ; it was evident that the comparison of EuESCs and EcESCs was masked by the large interpatient variation. Hence, we focused on the comparisons of oxygenation conditions and carried out separate normalization for EuESCs and EcESCs (i.e., data for EcESCs incubated in normoxia was also set to 100% for each individual PK; Figure 1B ). The obtained PK activity profiles revealed subtle yet statistically significant decrease in PKAc and Akt activity for EcESCs incubated in hypoxia vs normoxia (P < 0.05 and P ≤ 0.01, respectively). The trends were similar in the case of EuESCs incubated in hypoxia vs normoxia, yet statistical significance could not be achieved. No clear trends were observed in the activity of CK2, whereas in the case of ROCK, a slight elevation of kinase activity under hypoxic conditions was evident, which was however not statistically significant.

As the PK assay revealed only minor changes in the examined PK pathways which could not be explained solely by the effect of hypoxia, we decided to assess a wider picture of putative differences between the EuESC and EcESC signaling pathways under different oxygenation conditions by exploring the transcriptome.

To explore how hypoxia affects the gene expression of SCs isolated from paired EuE and EcE of women with endometriosis, we sequenced the entire transcriptome of the cell culture samples exposed to both normoxia and hypoxia (EuESCs, N = 5 and EcESCs, N = 5, in total 20 samples). By applying the gene filtering procedure as mentioned in the methods section, 17,887 genes were retained for subsequent analyses. The principal components analysis (PCA, Supplementary Figure 1 ) showed no specific grouping of samples based on the correlation between the data points representing each sample in four study groups. The Pearson correlation analysis revealed high intra-group variability in EuESCs exposed to hypoxia and EcESCs exposed to normoxia ( Supplementary Figure 2 ). However, there was a strong similarity at the transcriptome level between all the samples in the study groups based on the correlation coefficient values ranging from 0.85 to 1.

The comparisons of EuESCs and EcESCs exposed to hypoxia and/or normoxia yielded different numbers of differentially expressed genes (DEGs) (the comparison groups and numbers of DEGs are shown in Table 2 , statistically significant DEGs presented in Supplementary Table 1 ). The gene ontology and pathway enrichment analyses of DEGs in comparison groups identified statistically significantly enriched biological processes and pathways presented in Figure 2 and Supplementary Table 2 .

In EuESCs, we identified 116 DEGs upregulated under hypoxic conditions (log₂FC range from 0.5 to 2.9). Gene set enrichment analysis showed glycolysis/gluconeogenesis, biosynthesis of amino acids and carbon metabolism pathways among the top enriched pathways, which correlated with HIF-1 signaling pathway via shared DEGs ( Supplementary Figure 3A ). Among the top biological processes, the analysis revealed supramolecular fiber organization, NADH regeneration, glucose catabolic process to pyruvate, and canonical glycolysis ( Supplementary Figure 3B ). The expression profiles of EcESCs exposed to hypoxia vs normoxia were similar and the analysis identified a single statistically significant DEG, TGFBI (log₂FC 2.0).

In the comparison group of both EcESCs and EuESCs exposed to hypoxia, we identified 20 DEGs, of which 15 were downregulated (log₂FC range from –8.8 to –0.8) and 5 upregulated (log₂FC range from 0.7 to 3.4). The comparison of EcESCs and EuESCs both exposed to normoxia revealed 163 DEGs, 142 of which were upregulated (log₂FC range from 0.5 to 6.9) and 21 downregulated (log₂FC range from –0.5 to –2.3) in ectopic samples. Among the top enriched KEGG pathways in the upregulated DEG group, we identified focal adhesion, PI3K-Akt signaling pathway, and extracellular matrix (ECM)-receptor interaction ( Supplementary Figure 3C ). The anatomical structure morphogenesis, cell adhesion, cell migration, blood vessel morphogenesis, and blood vessel development were among the top enriched biological processes ( Supplementary Figure 3D ).

Given that EcE resides in vivo at more hypoxic conditions than the EuE, we hypothesized that in the case of EcESCs, incubation in a hypoxic environment closely resembles the physiological conditions, while in the case of EuESCs, hypoxia induces reductive stress. We therefore proceeded with a comparison of cell transcriptomes mirroring the native-like states, i.e., EcESCs exposed to hypoxia compared to EuESCs exposed to normoxia. We found 569 DEGs, with 375 upregulated genes (log₂FC range from 0.5 to 7.7) and 194 downregulated genes (log₂FC range from –0.5 to –8.7). The analysis of upregulated genes showed sugars and amino acids metabolism pathways, glycolysis, ECM-receptor interaction, focal adhesion, and PI3K-Akt signaling pathway and wound healing, blood vessel development, actin cytoskeleton development processes among the top enriched pathways and biological processes ( Supplementary Figures 3E, F ). In the downregulated DEGs group, RAS signaling and regulation of actin cytoskeleton were among the top KEGG pathways, and the regulation of cell differentiation and regulation developmental process were among the top biological processes ( Supplementary Figures 3G, H ).

The Venn diagram showed that the focal adhesion and the regulation of actin cytoskeleton KEGG pathways were shared among the three groups: EuESCs exposed to hypoxia vs normoxia, EcESCs vs EuESCs exposed to normoxia, and EcESCs exposed to hypoxia vs EuESCs exposed to normoxia ( Supplementary Figure 4 ). The analysis of common DEGs in the 5 comparison groups indicated that TGFBI was upregulated in all the comparisons ( Table 2 ; Figure 3 ).

According to previous studies, TGFBI is a secreted protein (44, 45). To establish whether the trends observed on the level of transcriptome are also valid at the protein level, we quantified the amount of secreted TGFBI in the spent cell culture media using ELISA. To compensate for the differences in the number of TGFBI-secreting cells following the 48-h incubation in different oxygenation conditions, the ELISA data was normalized to the total protein content in lysates of cells from which the media aliquots were collected. Furthermore, to compensate for the interpatient variation in TGFBI secretion, for each patient, the relative TGFBI protein content was additionally normalized according to that established for the EuESCs in normoxia (set to 100%).

The results (N = 5 in each group) are summarized in Figure 4A and Supplementary Figure 5 . Expectedly, relative TGFBI secretion was elevated in the medium collected from EuESCs incubated under hypoxic vs normoxic conditions (P < 0.01), and in the medium collected from EcESCs in both conditions relative to EuESCs (P < 0.001). In the case of EcESCs, no significant difference in TGFBI secretion was observed dependent on the oxygenation conditions, yet the trend was similar to that observed in transcriptomics with higher TGFBI levels in EcESCs exposed to hypoxia. Furthermore, in the case of samples collected from EcESCs in hypoxia, the TGFBI levels were closer to the upper limit of the calibration curve, thus potentially complicating the comparison.

As the isolated SCs represent a simplified model system, we proceeded to investigate the TGFBI in the context of endometriosis in a physiologically more relevant system represented by the full tissue samples. The ESCs used in our experiments were cultured in the absence of hormones to achieve a sufficiently numerous and homogenous cell population. As a result, the cells likely lost the characteristics of the menstrual cycle phase, at which the patient samples were initially collected (46). Therefore, we decided to investigate whether TGFBI mRNA expression varies depending on the menstrual cycle phase. For that, we determined TGFBI mRNA level in 45 EuE from women with endometriosis (11 in the proliferative phase and 34 in the secretory phase), and 24 EuE from women without endometriosis (controls, 5 in the proliferative phase and 19 in secretory phase) using qRT-PCR.

We found TGFBI mRNA expression in EuE from both women with endometriosis and controls to be menstrual cycle phase-dependent, with a higher expression level in the proliferative phase compared to the secretory phase (FC = 4.5, P < 0.0001 and FC = 4.5, P < 0.001, respectively, Figure 4B ). Having confirmed that TGFBI mRNA levels are dependent on the menstrual cycle phase, we proceeded to the exploration of TGFBI protein expression and localization in tissue samples from women with endometriosis obtained at proliferative or secretory phases.

To determine which cell types express TGFBI protein in endometrium and lesions, IHC was performed on paired samples of EuE and EcE from seven women with endometriosis, and additionally 17 EuE from women with endometriosis (samples collected during the proliferative or secretory phases of the menstrual cycle).

We observed a moderate and diffuse TGFBI protein expression in the stroma of paired EuE during the proliferative phase of the menstrual cycle ( Figure 4C ; Supplementary Figure 6 ). In the secretory phase, TGFBI protein staining in stroma was sparse and focal, localized around the vessels ( Figure 4D ; Supplementary Figure 6 ). The tissue sections of paired EcE contained fibrotic areas, and TGFBI protein expression was diffuse in stroma and around the vessels ( Figures 4E, F ; Supplementary Figure 6 ). In the secretory phase, TGFBI protein expression was stronger in EcE compared to EuE, which is similar to our observations in the SC culture experiments. Additionally, we analyzed TGFBI protein signal intensity in 17 EuE from women with endometriosis (6 in proliferative phase and 11 in secretory phase). TGFBI protein expression was significantly higher in the proliferative phase compared to the secretory phase (P < 0.05, Supplementary Figures 7 - 9 ). Due to the fibrotic nature of EcE, we did not compare TGFBI signal intensity between 6 paired EuE and EcE samples in the secretory phase to avoid potential artifacts.