Endometrial organoids derived from Mayer–Rokitansky–Küster–Hauser syndrome patients provide insights into disease-causing pathways

Mayer–Rokitansky–Küster–Hauser (MRKH) syndrome (OMIM: 277000) is a rare malformation, characterized by the partial or complete absence of the uterus and the upper two-thirds of the vagina, due to a still unknown defect in embryonic development (Oppelt et al., 2006). It affects 1 in 4500 women, making it the second most common reason for primary amenorrhea (Aittomäki et al., 2001; Herlin et al., 2016; Timmreck and Reindollar, 2003). Typically, MRKH patients have a normal female karyotype (46, XX) and regular development of secondary sexual characteristics, as their ovaries are functional. Although MRKH syndrome can occur as an isolated genital malformation (MRKH Type I), it is often associated with additional renal and/or skeletal, and to a lesser extent, with auditory, cardiac and other defects (MRKH Type II) (Oppelt et al., 2012; Rall et al., 2015b). In both cases, patients with MRKH frequently have one or two uterine remnants, consisting of myometrium and, although less often, also endometrium (Rall et al., 2013).

The etiology of the syndrome remains largely enigmatic, yet the spectrum of malformations encountered in MRKH patients suggests that the disease originates from a developmental defect during embryogenesis. Moreover, cases of familial clustering have implied a genetic component in the etiology (Herlin et al., 2014; Williams et al., 2017; Nik-Zainal et al., 2011). Mutations in WNT4, which is essential for the complete formation of the Müllerian ducts (MDs) (Vainio et al., 1999), have been reported earlier as a possible cause for MRKH in a small number of patients (Biason-Lauber et al., 2004; Philibert et al., 2008, 2011). Since WNT4 is also necessary to prevent formation of Leydig cells in women, patients with mutated WNT4 also present with clinical hyperandrogenism, rendering it a slightly different clinical entity (Biason-Lauber et al., 2004, 2007; Philibert et al., 2008). Furthermore, in some patients, mutations and possibly harmful variants have been found in developmental genes like WNT9B (Wang et al., 2014; Waschk et al., 2016), LHX1 (Ledig et al., 2011, 2012; Sandbacka et al., 2013) or TBX6 (Sandbacka et al., 2013; Tewes et al., 2015, 2019). Of high interest are specifically LHX1, WNT4 and WNT9B due to their previously reported role in the formation of the MDs from the coelomic epithelium in gestational week 6 (Mullen and Behringer, 2014). After the two MDs are formed, they start growing caudally along the Wolffian ducts (WDs). By week 8, the terminal ends of the MDs begin to fuse, forming the uterovaginal duct. In males, the MDs start to regress after week 10 under the influence of anti-Müllerian hormone (AMH) and WNT7A. In females, however, they differentiate into the fallopian tubes, uterus, cervix and vagina, specifically regulated by the coordinated action of transcription factors and signaling molecules such as homeodomain transcription factors (e.g. Hox genes) and members of the WNT family (Robboy et al., 2017; Roly et al., 2018). This concerted interplay between transcription factors, hormones and growth signals during embryogenesis leads to a variety of highly timed and spatial gene expression changes (Roly et al., 2020). The widespread lack of a clear genetic link in MRKH syndrome suggests that the answer might be found at the transcriptional level. Therefore, the focus in recent years has shifted towards molecular characterizations of primary diseased tissue, identifying a plethora of new candidate genes and regulatory networks that might drive or contribute to the pathology (Hentrich et al., 2020; Rall et al., 2015a). In addition, hormonal stimulation of primary endometrial stroma cultures derived from MRKH rudimentary tissue showed a reduced transcriptional response compared to healthy controls (Brucker et al., 2017), indicating that dysfunctional hormone receptors play a role in the pathophysiology of MRKH (Rall et al., 2013). However, as other estrogen- and progesterone-dependent tissues like the breast develop normally in MRKH patients, any possible defect would have to be limited to the genital tract. This stresses the importance of studying the rudimentary tissue directly. Whereas stromal cultures have already been investigated (Brucker et al., 2017), the absence of a functional model for the glandular epithelium of MRKH endometrium to analyze the characteristics and capabilities of affected cell types vastly limited the pathophysiological understanding of this disease. Recent years have seen the establishment of patient-derived organoid models from healthy and diseased endometrium (Turco et al., 2017; Boretto et al., 2017, 2019). Organoids are self-renewing, three-dimensional (3D) models that mimic key properties and characteristics of the original in vivo tissue, which greatly facilitates research into complex interactions of the involved cell types.

Here, we show for the first time the successful establishment of organoid models from endometrium found in the rudimentary tissue of MRKH patients. The established organoids showed high similarity to healthy endometrial organoids, were hormone responsive and could be cultured long term. Yet, gene expression analyses by RNA sequencing showed distinctive differences between diseased and healthy organoids, emphasizing their potential as a powerful tool to investigate the etiology of MRKH syndrome. We identified several important developmental transcription factors to be differentially expressed between healthy and MRKH organoids. This highlights the importance and value of these models for studying and understanding the pathogenesis. The establishment and long-term growth of epithelial tissue from MRKH rudiments provides a powerful addition to the toolbox for studying the disease in a controlled and standardized environment, and these endometrial organoids promise to pave new avenues towards a better understanding and possible treatments of the disease.

Human endometrial organoid cultures have recently been established by several research groups (Turco et al., 2017; Boretto et al., 2019, 2017). Typically, glandular structures obtained after processing endometrial biopsies are embedded in an extracellular matrix (ECM) component such as basement membrane extract (BME) or Matrigel and cultured in the presence of a defined cocktail of growth factors (Table S2). Since endometrium can be found in MRKH uterine rudiments (Rall et al., 2013), we screened a cohort of MRKH patients that underwent a laparoscopically assisted creation of a neovagina for the presence of endometrium in uterine rudiments. A total of 48 patients [35 MRKH Type I (73%), 13 MRKH Type II (27%)] were screened, of which 37 [32/35 Type I (91%), 5/13 Type II (38%)] had uterine rudiments present (Fig. 1A). We were able to macroscopically detect endometrium (Fig. 1B,C) in 12 patients [12/32 Type I (33%), 0/5 Type II (0%)]. Endometrium of four rudiments was also subjected to histological processing (Fig. 1B, asterisk) and immunohistochemistry for the transcription factor PAX8 (highly expressed in endometrial epithelial cells), estrogen receptor alpha (herein referred to as ER; encoded by ESR1), progesterone receptor (PR; encoded by PGR) and the proliferation marker Ki67 (also known as MKI67) (Fig. 1D) to verify the endometrium entity of the tissue. PAX8 staining confirmed that the endometrial epithelium of MRKH patients had all the morphological features of a normal endometrium, showing tubular and frequently branching glands with a single layer of columnar epithelium (Fig. 1D), as it has been shown recently in an extensive histology study of MRKH rudiments (Rall et al., 2013). Whereas ER and PR were ubiquitously expressed in the epithelial as well as stromal cells of the endometrium, the proliferative capacity, measured by Ki67 expression, was almost absent in the epithelium and stroma of MRKH patients (Fig. 1D), as reported previously (Rall et al., 2013). Although the proliferative capacity of the initial cell population was low, we successfully established 14 organoid models (ten patients had endometrium present unilaterally and two patients bilaterally, which accounts for the two additional organoid models) (Table S1). The success rate in establishing organoid models from patients with a macroscopically detected endometrium was 100% (14/14). The fact that the endometrium size of MRKH patients ranged from a few millimeters in diameter (e.g. MRKH #03, Fig. 1D) to even under a millimeter (e.g. MRKH #04, Fig. 1D) limited the amount of starting material for culture setup, but did not hinder the successful establishment of an organoid model (compare MRKH #03 and #04, Fig. 1E). Within 3–6 days, cystic organoid structures were visible, and all 14 established MRKH organoid models as well as four control models (Fig. 1E; Fig. S1A,B) were successfully cultured long term for more than 15 passages (>6 months of culture) without showing signs of a decrease in proliferation (Fig. S1C). Notably, the organoid models could easily be cryopreserved and stored as a live biobank (see Materials and Methods).

The endometrial organoids of MRKH patients and healthy controls were expanded and passaged at ratios of 1:2 or 1:3 every 10–14 days and characterized by immunohistochemistry and immunofluorescence (Fig. 2A). Similar to the primary tissue, the endometrial organoids showed high and ubiquitous expression of PAX8 and ER, whereas the absence of steroid hormones (e.g. estradiol) in the culture medium explains the absence of PR expression in all models. The addition of estradiol to the culture medium led to the restoration of PR expression in both control and MRKH organoids (Fig. S1D). In sharp contrast to the observations in MRKH endometrial tissue, the organoids of MRKH patients expressed high levels of the proliferation marker Ki67, comparable to healthy controls. Markers of glandular epithelium [cytokeratin, epithelial cell adhesion molecule (EpCAM) and E-Cadherin] were ubiquitously expressed in all organoid models, and Perlecan staining at the basolateral membrane throughout the entirety of the endometrial organoids showed that epithelial polarity remained intact (Fig. 2B). Interestingly, phenotypic and morphological characterizations of healthy and diseased organoid models revealed no obvious differences.

After investigating organoid morphology, we next interrogated their transcriptome. Using RNA sequencing, we profiled seven organoid models from MRKH patients and four from healthy controls that were grown in expansion medium (ExM), treated with beta-estradiol (E2) or the combination of beta-estradiol and progesterone (E2+P4). Based on these six experimental groups, differential gene expression was determined according to disease status and treatment with cut-offs of Benjamini–Hochberg-adjusted P-value (P_adj)≤0.05 and |log₂ fold-change (FC)|≥0.5 (Fig. 3A; Table S6).

In the principal component analysis (PCA) of expression profiles, the samples partitioned clearly along principal components 1 (PC1) and 2 (PC2). The grouping along PC1 into ExM-treated samples to the right and hormone-treated samples to the left indicated the strongest influence on expression differences through treatment (Fig. S2). Accordingly, the expression effects of E2 and E2+P4 treatment led to about five times more differentially expressed genes (DEGs) compared to disease status comparisons, which yielded 492 DEGs for ExM treatment (Fig. 3A). The PCA further suggested great homogeneity in organoid responses with respect to treatment as reflected by the sample partitioning along PC2 according to patient origin (Fig. S2). Moreover, computational estimation of the underlying cell-type composition through single-cell reference data from uterus (Wang et al., 2020) and endometrial epithelial organoids (Fitzgerald et al., 2019) showed very consistent signatures across all samples with highest expression for epithelial cell types, agreeing with the model and in line with immunostainings for epithelial markers (Fig. 2B; Figs S3 and S4).

Focusing on differential expression under ExM first, we identified 492 DEGs with 365 upregulated and 127 downregulated genes (Fig. 3B). Among them were 12 homeobox proteins (upregulated in MRKH: LHX1, HOXD8, ONECUT3, LBX2, HOXB4, SATB1 and HOXB6; downregulated in MRKH: EMX2, ZHX3, IRX3, NKX6-2 and IRX5). Of those, LHX1, EMX2 and the Hox genes have previously been associated with MRKH syndrome either by sequencing of patient blood or functional in vivo studies with animal models (Ledig et al., 2012; Masse et al., 2009; Miyamoto et al., 1997). Moreover, supporting the role of the 492 DEGs in MD development, a large proportion of these genes are specifically activated during duct morphogenesis in chicken (Fig. S5) (Roly et al., 2020).

Applying enrichment analyses to identify affected pathways and cellular processes as well as transcriptional regulators potentially driving the differential expression, we identified ‘plasma membrane’ and ‘cell periphery’ as the most significant Gene Ontology (GO) terms and ‘ECM-receptor interaction’ with respect to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Fig. 3C). A binding site analysis suggested highly significant overrepresentation of the SP4 motif among DEG promotors (Fig. 3D). Interestingly, expression of SP4 was decreased in MRKH organoids (Fig. 3E), agreeing with previous observations of a smaller uterus in Sp4-knockout mice (Göllner et al., 2001). Complementary analyses that utilized chromatin immunoprecipitation-sequencing (ChIP-seq) data and thereby also accounted for indirect binding events, as well as transcription factors with less clear motifs (Puente-Santamaria et al., 2019), suggested the DEG set to be enriched for EZH2 and SUZ12 targets (Fig. 3F). Since EZH2 was also identified in our previous study of MRKH endometrial tissue (Hentrich et al., 2020), we saw accumulating evidence that further suggested epigenetic investigations into the origins of the syndrome.

Similarity between MRKH organoids and the previously analyzed endometrial tissue also existed with respect to DEGs. Comparing the set of 492 DEGs identified in MRKH organoids to the set of 2121 DEGs previously reported for MRKH endometrial tissue (Hentrich et al., 2020) led to 86 shared genes, of which 51 were dysregulated similarily in both organoid and tissue (Fig. 4A). Among them were the GATA5 transcription factor as well as the HOXB4 and HOXB6 homeodomain proteins, which were upregulated in both MRKH tissue and organoid (Fig. 4B).

The relatively small overlap of DEGs between diseased tissue and organoid might partially be attributable to differences in tissue composition. Of the 2121 DEGs originating from the mixture of cell types (stromal, epithelial, endothelial and blood) in primary patient tissue (Hentrich et al., 2020), about one-third showed virtually no expression in organoids based on epithelial cells alone. Nevertheless, the overlapping DEGs that do exist suggest that the organoids capture parts of the pathology in a highly homogeneous and reproducible way.

Next, we sought to investigate transcriptional changes upon treatment with steroid hormones in endometrial organoids of patients and controls, as the human endometrium undergoes substantial remodeling mainly controlled by the steroid hormones E2 and P4 during the menstrual cycle (Aghajanova et al., 2008; Roy and Matzuk, 2011). In order to investigate the hormone response of disease versus control organoids, RNA-sequencing data of E2- and E2+P4-treated organoids were analyzed for gene expression changes.

Under E2 treatment, 2321 DEGs were identified in control and 2290 DEGs in MRKH organoids (Fig. 5A), of which about three-quarters overlapped (Fig. 5B). Furthermore, 594 of the DEGs have been described previously for endometrial epithelial organoids treated with E2 (Fitzgerald et al., 2019). As already indicated by the PCA (Fig. S2), the expression changes were highly homogenous between samples (Fig. S6A) and encompassed genes such as FOXJ1, which activates essential genes for motile cilia formation and function, as well as DYDC2, a marker for ciliated cells (Fig. S6B). Consistently, the most significant GO terms for these DEGs were ‘cilium organization’ and ‘cilium assembly’ (Fig. 5C); a fact that was also observed by immunostaining of hormone-treated organoids with the cilia marker acetylated alpha-tubulin (Fig. 5D). In addition, the expression of genes attributed to ciliated cells strongly increased upon treatment with E2 in both MRKH and control organoids (Figs S3 and S4).

Intriguingly, despite largely similar gene expression changes upon E2 treatment between MRKH and control organoids, a few genes were specific to disease condition (Fig. 5E). Among them were genes such as AMY1 (also known as AMY1A) with upregulation upon E2 treatment specifically in MRKH organoids (Fig. 5F). Consistently and previously unnoticed, AMY1 also showed increased expression in MRKH endometrial tissue (Fig. S6C). Hence, amylases – typically not expressed in endometrial tissue – might have been spuriously upregulated in MRKH patients, leading to tissue breakdown and degeneration.

As the next step, we studied the response of MRKH and control organoids to the combination of E2 and P4. As mentioned above, the DEG count was very similar and symmetric between E2 alone and the combination of E2 and P4 (Fig. S7). In fact, for both control (Fig. S7A) and MRKH organoids (Fig. S7B), gene expression changes under the hormone treatments were almost identical as reflected by the nearly perfect standard diagonal (Fig. S7C,D), and very few genes showed deviations (Fig. S7E), indicating that the addition of P4 to E2 led to few additional transcriptome changes.

Finally, we validated disease- and patient-specific gene expression changes of the selected genes using quantitative PCR in the RNA-sequencing cohort as well as in an independent cohort (Fig. 6A). We focused on the DEGs that are also highly expressed: LHX1 (associated in Müllerian agenesis; Huang et al., 2014), HOXD8 (homeobox gene that is highly expressed during the development of the chicken MD; Roly et al., 2020), FAM3B (a recently identified FGFR ligand implicated in posterior development; Zhang et al., 2021), NDN, androgen receptor (AR) and GATA5 (expressed during the MD development and associated with abnormalities of the genitourinary tract of female mice; Roly et al., 2020; Molkentin et al., 2000). Quantitative PCR assays clearly validated the RNA-sequencing results (Fig. 6B). Additionally, the expression changes found in our MRKH sequencing cohort were also found in the three MRKH models that were selected as an independent cohort (MRKH-IC) (Fig. 6B).

During the last decade, a plethora of research projects, mainly in the field of genetics, have aimed at deciphering the cause of the multi-faceted MRKH syndrome (Herlin et al., 2020; Fontana et al., 2017). Reports of familial cases alongside seemingly sporadic appearances have fueled hope to find a genetic cause of the disease. Yet, a common genetic determinant has remained elusive thus far. Sequencing studies are mainly performed on blood samples and the focus lays primarily on mutations in the patients' germline. That, however, negates the possibility that the origin of MRKH syndrome lies in the developing tissue. MD development (Santana Gonzalez et al., 2021) is, as are all developmental processes, a highly complex sequence of events that involves activation and silencing of signaling pathways, expression of transcription factors (e.g. homeobox proteins) and cascades of growth factors, all of which must be precisely coordinated in space and time to arrive at the desired outcome (Mammoto and Ingber, 2010). Disturbances in these pathways either by lacking factors or their deregulated expressions may, hence, lead to developmental defects, as seen in MRKH patients with uterus aplasia alone (MRKH Type I) or with associated malformations (MRKH Type II).

To find potential disease-causing alterations, we strongly advocate for focusing on diseased tissue, rather than sequencing patient blood. So far, only one study reports the successful culturing of endometrial stromal cells isolated from MRKH rudimentary uterine horns (Brucker et al., 2017). Here, we were able to isolate the endometrium from rudimentary uterine horns and establish long-term proliferating epithelial organoid cultures that share high similarity with healthy controls and are fully hormone responsive (Figs 2 and 5). This is in stark contrast to the previously isolated endometrial stromal cells that showed impaired decidualization upon steroid hormone exposure (Brucker et al., 2017). We hypothesize that these results point to a pivotal interplay between stromal and epithelial cells. In vivo, estrogen-stimulated stromal cells secrete growth factors (e.g. EGF and FGF-10) that lead to the proliferation of epithelial cells (Chung et al., 2015). Hence, impairment of estrogen receptivity in the stromal compartment might, in turn, lead to a reduction of these factors in MRKH patients and result in the reduced proliferation capacity seen in vivo, but not in the organoid model due to the supplementation of these factors in the growth medium (Table S2). Successful establishment and proliferation of endometrial cultures show that these cells are intrinsically able to proliferate once they are provided with external growth signals. These results suggest an essential crosstalk and interplay between stromal and epithelial cells in the pathogenesis of MRKH. The organoid model established here lends itself to further explore this theory and to continue to explore the causal factors of the disease.

During development of the female reproductive tract, a cascade of tightly regulated processes in the developing embryo need to take place. This is largely governed by tightly controlled transcription factor programs, including homeobox-domain proteins. Clinically, MRKH represents a heterogeneous disease and is generally categorized into Type I (isolated finding) or Type II (accompanied by abnormalities of additional organ systems, primarily including the kidneys and the skeleton) (Herlin et al., 2020). Interestingly, in our patient cohort, about 90% of MRKH Type I cases had rudimentary horns, in contrast to Type II cases with only 38% (Table S1). The multifaceted appearance of the uterine horns (absent, one versus two, with endometrium or without and size variation) might point to defects at different phases of uterine development.

The uterus, both fallopian tubes, cervix and the upper third of the vagina develop from the MDs. This development can be divided into three phases (Santana Gonzalez et al., 2021). The first phase is the specification phase when Müllerian precursor cells are specified in the surface epithelium at the cranial pole of the mesonephros (embryonic kidney). During the second phase, termed invagination, these cells start to become mesoepithelial and move between the mesenchyme in a caudal direction. In the final phase, the elongation phase, invaginating cells from the MD epithelium contact the WD and the duct elongates. At this point, the most caudal part of the MD fuses into a single tube with a separating midline epithelial septum, forming the midline uterovaginal canal between the two adjacent MDs, which later disappears. The degree of midline fusion is vastly different between species. In mice, most of the MDs remain unfused and later give rise to large bilateral uterine horns instead of a single uterus in humans (Spencer et al., 2005). At this stage of development, the embryo is sexually indifferent and presents with both WDs and MDs (Kobayashi and Behringer, 2003). Under the influence of several signaling cascades and factors (e.g. AMH), the MDs regress in males, whereas lack of these hormones in females leads to regression of the WDs (Kobayashi and Behringer, 2003). After female sex determination, expression of Hox genes and retinoic acid (RA) signaling drive segmentation of the proximal and distal MD (Santana Gonzalez et al., 2021). Each specified Müllerian segment will give rise to a different part of the female reproductive tract, with the most anterior differentiating into the oviducts, and the most posterior into the upper vagina. We hypothesize that there is not a single defect in this finely tuned developmental process that causes MRKH syndrome, and that the multifaceted clinical phenotype can be attributed to the fact that different pathways were altered in different subsets of patients. MRKH patients without rudiments and patients with two developed rudiments containing endometrium were very likely to have different alterations in their genome and/or epigenome. The latter seemed to have had a caudal MD fusion defect whereas the former phenotype probably arose from a defect during the elongation phase. The fact that MRKH patients have fallopian tubes (Herlin et al., 2020), both emerging from MDs, provides evidence that the initial steps of MD development are still unchanged. Another possibility might be that the regression of the WDs (driven by apoptosis and degradation of tissue, which in turn are driven by metalloproteinases; Page-McCaw et al., 2007) in MRKH patients malfunctioned and led to partial regression of the adjacent MDs, which in turn led to a fusion defect in the uterovaginal canal. Hence, the structure regressed and could not establish a contact with the urogenital sinus that was forming the lower part of the vagina. Investigation of endometrial tissue of patients with other Müllerian anomalies such as bicornuate uterus, uterus didelphys (defect in fusion) or a uterus septum (defect in septal resorption) might provide important insights in the causality of MRKH as these likely represent similar, but milder defects, compared to MRKH patients.

The ECM has been shown to play a pivotal role in various biological processes and is essential for the function and integrity of animal tissues (Lu et al., 2011). During development, the ECM is constantly being remodeled by degradation and reassembly processes. The female reproductive tract development, as outlined before, is marked by numerous such events of ECM remodeling. Interestingly, our RNA-sequencing results show an increase in the expression of genes involved in specific cellular processes and pathways related to the cellular compartments ‘plasma membrane’ and ‘cell periphery’, as well as ‘ECM-receptor interaction’, which were found to be highly enriched (Fig. 3C; Table S5). Although the GO terms ‘plasma membrane’ and ‘cell periphery’ are primarily only an indication of the location of the deregulated proteins, they suggest that the majority of DEGs in endometrial epithelial cells in MRKH organoids as compared to healthy controls are implicated in the interaction of glandular epithelial cells with the ECM and/or the stromal compartment of endometrium in the uterine rudiments. Moreover, when grown in culture, epithelial cells from MRKH patients have the potential to grow, respond to hormones and, in fact, behave very similarly to healthy organoids. This hints at the importance of the communication between epithelial and their surrounding cells for the normal development of the uterus, and highlights the importance of further investigating these interactions. Complex interactions between different cell types of the MRKH uterine rudiments have thus far not been implemented in the research community. Isolating stromal cells, as previously reported (Brucker et al., 2017), and performing co-culture experiments with endometrial organoid models might be a promising way of tackling this highly understudied research field. Furthermore, investigations using mouse models might be a promising tool to verify whether altered gene expression or knockouts of the candidate genes identified here might lead to genital malformations such as those seen for LHX1 and other candidate genes (Herlin et al., 2020).

Uterine development occurs during the first trimester of development. At the time of surgery and tissue sampling, our patients have passed this point on average by 20 years (with one patient already being 41 years old) (Table S1). It is remarkable that we were able to establish fully functional endometrial organoids from this tissue, and even more remarkable that we were able to find specific expression differences in 492 DEGs including already suspected candidate genes such as LHX1 and Hox genes. Surprisingly, these were mainly found to be expressed in MRKH organoids and absent in healthy controls. We speculate that this explains why any single candidate gene may not necessarily be mutated in all patients. Mutations that cause a complete loss of function might have very severe developmental defects leading to embryonic lethality and, hence, cannot be discovered. For example, homozygous Lhx1 knockout in mice leads to only a few neonates that are born with an MRKH-like phenotype, but also lack anterior head structures (Kobayashi et al., 2004). A deregulation (but not necessarily a knockout) of this gene, specifically during uterine development or in specific regions of the developing female reproductive tract, might however cause MD anomalies while leaving other regions unharmed. The fact that a large portion of the MRKH candidate gene identification data is curated from knockout studies in animal models (Masse et al., 2009) might be the reason why many of our identified genes have not yet been associated with the disease. The transcription factor GATA5, for example, is highly expressed in MRKH endometrial organoids as well as in tissue (Hentrich et al., 2020) compared to healthy controls, but has never been linked with MD anomalies, even though it is highly and dynamically expressed in the developing tissue, as revealed by a large-scale RNA-sequencing study in chicken (Roly et al., 2020). In the future, shifting the focus from knockout studies towards interrogation of dysregulated gene expression using CRISPR-A screens (Kampmann, 2018) might help to better understand the processes during uterine development. Several of our identified DEGs could serve as a starting point to select interesting candidates.

MRKH syndrome is the most common cause of uterine aplasia at a frequency estimated to be 1 in 4500 female newborns. Since the early 2000s, treatment of MRKH patients mainly focuses on the co-existing vaginal aplasia by creation of a neovagina to give patients the possibility of sexual intercourse (Rall et al., 2014). This vastly improves the quality of life of affected women (Weijenborg et al., 2019). However, this treatment is not adequate to provide a solution for the missing reproductive ability in MRKH patients. Besides adoption and surrogacy, uterus transplantation has emerged as a highly promising method to overcome this shortcoming (Brannstrom et al., 2018). The two latter methods come with a plethora of ethical difficulties, and researchers and clinicians are working urgently to find new ways of providing larger groups of patients with a cure. Organ-on-a-chip-based methods as well as using animals or deceased patient uterus scaffolds represent opportunities for uteri reconstruction using patient-specific cells to overcome limitations of transplantation approaches (Bergmann et al., 2021). The development of these techniques in combination with an in-depth understanding of uterine development might mark a milestone for achieving a long-term treatment solution to some of the most severe uterine pathologies such as MRKH syndrome. To use patient-derived cells for uterus reconstruction, we need to get a deeper understanding of their functionality and their interactions. Hence, endometrial organoids from MRKH patients and their interaction with stromal cells need to be investigated further to achieve this goal.

Our study shows that MRKH endometrial organoids are hormone responsive and show a high similarity to healthy endometrial epithelial cells. Future research in this emerging line of research, however, should reveal whether this is sufficient for successful embryo implantation in a to-be-developed de novo uterus from MRKH patient cells (Alzamil et al., 2021; Bergmann et al., 2021).

All tissue biopsies were obtained from patients after informed written consent. The study was conducted according to the guidelines of the Declaration of Helsinki and was approved by the Ethics Committee of the Eberhard Karl University of Tübingen (Ethical approval 205/2014BO1 and 150/2018BO2) and is compliant with all relevant ethical regulations regarding research involving human participants. Rudimentary uterine horns were excised from 37 MRKH patients at the time of laparoscopically assisted creation of a neovagina and transported to the pathology lab in sterile containers. Uterine rudiments were sectioned perpendicular to the longest axis and macroscopically evaluated for presence of endometrial tissue. If present, endometrium was removed with as little attached myometrium as possible and further processed. For this study, endometrial tissue from 12 patients with MRKH syndrome was collected. Endometrial biopsy samples [obtained via Pipelle Endometrial Suction Curette (Medesign IC)] of four premenopausal patients served as controls. For full patient characteristics, see Table S1.

Tissue samples were minced into small pieces (1–3 mm³) using a scalpel and dissociated with collagenase/dispase (1 mg/ml; COLLDISP-RO, Roche) in the presence of Rock inhibitor (Y-27632, 10 μM; M1817, Abmole Bioscience) for 1 h at 37°C on a shaking table. The digestion was attenuated by addition of medium [advanced Dulbecco's modified Eagle medium (DMEM)/F12 (12634010, Gibco) without serum] and centrifuged at 478 g for 10 min. The final pellet was resuspended in advanced DMEM/F12 supplemented with 1% Glutamax, 1% HEPES and 1% penicillin/streptomycin (all Gibco) and the desired amount of cell suspension mixed with BME (Type II, 3533-001-2, Trevigen) at a ratio of 65% BME to 35% cell suspension. 20 µl droplets were plated out on pre-warmed 48-well plates and placed upside-down in an incubator (37°C, 5% CO₂) for solidification. Afterwards, culture medium (Table S2) was added to each well and renewed every 3 days. Noggin conditioned medium from HEK293T-Noggin-Fc-cells (kindly provided by Hans Clevers, Hubrecht Institute, Utrecht, Netherlands) was produced as previously described (Farin et al., 2012). R-Spondin conditioned medium was produced with Cultrex^® HA-R-Spondin1-Fc 293T Cells (Trevigen) according to the distributor's protocol.

For passaging or cryopreservation, organoids were recovered by resuspending the BME drops in ice-cold advanced DMEM/F12 and transferred to 15 ml tubes. The organoid suspension was either mechanically or enzymatically [25% 1× TrypLE Express (Gibco), 75% advanced DMEM/F12] dispersed and then pelleted. For further culture, the pellet was reconstituted in advanced DMEM/F12 and mixed with BME at a ratio of 65% BME to 35% cell suspension and cultured as described above. For cryopreservation, the cell pellet was resuspended in Recovery™ Cell Culture Freezing Medium (Gibco), the solution transferred to cryo-vials and then cooled down in CoolCell™ LX Freezing Containers (Merck) in a −80°C freezer. The next day, vials were transferred for long-term storage to liquid nitrogen tanks.

Tissues and organoids were fixed in 4% paraformaldehyde followed by dehydration, paraffin embedding, sectioning (4 µm) and standard Haematoxylin and Eosin staining. Immunohistochemistry was performed on a Ventana Discovery automated immunostaining system (Ventana Medical Systems, Tucson, AZ, USA) using antibodies as specified in Table S3.

For characterization of organoids, formalin-fixed paraffin-embedded (FFPE) sections (4 µm) were subjected to heat-induced antigen retrieval and incubated with primary antibodies as specified in Table S3. Sections were incubated overnight in a humidified chamber at 4°C, and afterwards washed three times with PBS containing 0.1% Tween-20 and incubated for 1 h at room temperature with the respective secondary antibody (Table S3). Finally, the sections were again washed and mounted in ProLong Diamond Antifade mounting media containing DAPI (Thermo Fisher Scientific). Images were acquired with the EVOS M7000 imaging system and processed using ImageJ (NIH, USA).

Organoids were passaged as described above and allowed to grow for 4 days in standard culture medium (expansion medium; ExM). Afterwards, three different groups were established. Group 1 (ExM) served as an untreated control sample and received ExM for an additional 6 days, with the medium refreshed every other day. Group 2 (E2) was cultured with ExM containing 10 nM beta-estradiol (E8875, Sigma) for 6 days. Again, the medium was renewed every other day. Group 3 (E2+P4) first received ExM supplemented with 10 nM beta-estradiol for 2 days and afterwards ExM with 1 µM progesterone (P8783, Sigma) and 1 µM cAMP (1140, Tocris) in addition to 10 nM beta-estradiol for additional 4 days. After a total of 10 days, all organoids were harvested. For this, medium was removed from the wells and the BME domes incubated with 1× TrypLE Express at 37°C for 10 min. PBS was added to dilute the TrypLE Express and the suspension centrifuged at 478 g for 10 min. After discarding the supernatant, the pellet was resuspended in PBS and centrifuged again to remove leftover BME. After removing the supernatant, cell pellets were snap frozen in liquid nitrogen and then stored at −80°C.

Total RNA was isolated from organoid cultures using the RNeasy Mini Kit (74004, Qiagen). Simultaneous elimination of genomic DNA was achieved with on-column DNA digestion (79254, Qiagen) according to the manufacturer's protocol. The quality was assessed with an Agilent 2100 Bioanalyzer. Using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina and 100 ng of total RNA for each sequencing library, poly(A)-selected paired-end sequencing libraries (101 bp read length) were generated according to the manufacturer's instructions. All libraries were sequenced on an Illumina NovaSeq 6000 platform at a depth of around 40 million reads each. Library preparation and sequencing procedures were performed by the same individual, and a design aimed to minimize technical batch effects was chosen.

Read quality of RNA-seq data in fastq files was assessed using FastQC (v0.11.9) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to identify sequencing cycles with low average quality, adaptor contamination or repetitive sequences from PCR amplification. Reads were aligned using STAR (v2.7.6a) (Dobin et al., 2013) allowing gapped alignments to account for splicing against the Homo sapiens genome from GENCODE v35. Alignment quality was analyzed using SAMtools (v1.10) (Li et al., 2009). Normalized read counts for all genes were obtained using DESeq2 (v1.32.0) (Love et al., 2014). Transcripts covered with less than 50 reads (median of all samples) were excluded from the analysis, leaving 15,413 genes for determining differential expression. Cut-offs of |log₂FC|≥0.5 and P_adj≤0.05 were set to determine differentially expressed genes. Gene-level abundances were derived from DESeq2 as normalized read counts, and used for calculating the log₂-transformed expression changes underlying the expression heatmaps, for which ratios were computed against mean expression in control samples. The normalized counts provided by the ‘sizeFactors’ function in DESeq2 also went into calculating normalized reads per kilobase per million total reads (nRPKMs) as a measure of relative gene expression (Srinivasan et al., 2016).

g:Profiler2 (v0.2.0) was employed to identify overrepresented GO terms for differentially expressed genes (Raudvere et al., 2019). Transcription factor-binding site analyses were carried out in Pscan (v1.6) (Zambelli et al., 2009) on the Homo sapiens genome considering the −450 to +50 bp of the promoter region for motifs against the Jaspar 2020_NR database. TFEA.ChIP (v1.12.0) was employed with default parameters to determine transcription factor enrichment using the initial database version of ChIP-seq experiments (Puente-Santamaria et al., 2019). Cell-type-specific endometrial marker genes were taken from a study with single-cell profiles from uterine tissue (Wang et al., 2020), as well as endometrial epithelial organoids (Fitzgerald et al., 2019).

Total RNA was extracted from organoids using the RNeasy Mini Kit (74004, Qiagen), simultaneously eliminating genomic DNA with on-column DNA digestion (79254, Qiagen). Equal amounts of total RNA (1 μg) were reverse transcribed using the Maxima™ H Minus cDNA Synthesis Master Mix (M1661, Thermo Fisher Scientific). To analyze gene expression, 5 ng cDNA was subjected to real-time qPCR using PowerUp™ SYBR^® Green Mastermix (A25741, Thermo Fisher Scientific) on the QuantStudio 3 Real-Time PCR system (A28567, Thermo Fisher Scientific). Thermal cycling was performed with 3 min at 95°C, followed by 40 cycles of 95°C for 15 s and 60°C for 60 s. The specificity of the RT-qPCR products was assessed by melting curve analysis. Relative quantification was performed using the 2^−ΔΔCt method with SDHA and RPL13A as reference genes. Expression was normalized to the endometrial control group. All experiments were performed in duplicates. PCR primers are listed in Table S4.

We thank all patients who participated in the study. Special thanks go to Ingrid Teufel and Sabine Hofmeister for excellent technical support. We thank the team of pathologists of the Department of Anatomic Pathology, University of Tübingen and the NGS Competence Center Tübingen, foremost Nicolas Casadei and Olaf Rieß, for their support. We thank Nico Weber for his input in the initial analysis of the RNA-sequencing data. Special thanks go to Ana Maia dos Santos Leite for her critical input throughout the study and for her valuable comments on the manuscript.