Molecular Landscape of Pelvic Organ Prolapse Provides Insights into Disease Etiology

Pelvic organ prolapse (POP) represents a major health care burden in women, but its underlying pathophysiological mechanisms have not been elucidated. We first used a case-control design to perform an exome chip study in 526 women with POP and 960 control women to identify single nucleotide variants (SNVs) associated with the disease. We then integrated the functional interactions between the POP candidate proteins derived from the exome chip study and other POP candidate molecules into a molecular landscape. We found significant associations between POP and SNVs in 54 genes. The proteins encoded by 26 of these genes fit into the molecular landscape, together with 43 other POP candidate molecules. The POP landscape is located in and around epithelial cells and fibroblasts of the urogenital tract and harbors four interacting biological processes—epithelial-mesenchymal transition, immune response, modulation of the extracellular matrix, and fibroblast function—that are regulated by sex hormones and TGFB1. Our findings were corroborated by enrichment analyses of differential gene expression data from an independent POP cohort. Lastly, based on the landscape and using vaginal fibroblasts from women with POP, we predicted and showed that metformin alters gene expression in these fibroblasts in a beneficial direction. In conclusion, our integrated molecular landscape of POP provides insights into the biological processes underlying the disease and clues towards novel treatments.


Introduction
Pelvic organ prolapse (POP) represents a major health care burden in women, with a reported prevalence of up to 40% [1]. POP is not life-threatening, but its bothersome symptoms have a significant impact on quality of life. The cumulative life-time risk for a surgical intervention for POP is approximately 13% [2], and the total costs for POP surgery are substantial [3]. POP has a multifactorial etiology in which both hereditary and environmental factors play a role. Parity, vaginal delivery, birth weight, age, and high body mass index (BMI) have repeatedly been identified as environmental risk factors [4,5] and could all cause defects of pelvic floor tissues. Familial susceptibility to these defects seems likely, since women with POP more often have family members with POP compared to women without POP [6][7][8]. Evidence for a genetic background of POP has also emerged from twin studies, with approximately 40% of the POP liability being explained by genetic factors [9,10]. Candidate gene association studies have found a number of single nucleotide polymorphisms (SNPs) to be putatively associated with POP, including genes encoding steroid hormone receptors, collagens type I and III, and multiple matrix metalloproteinases. However, only one SNP, rs1800012 in the COL1A1 gene, has demonstrated a replicable association with POP [11]. Thus far, two genome-wide association studies (GWASs) of POP have shown significant associations with SNPs. The first GWAS identified five SNPs associated with POP that were nominally replicated [12]. In 2020, the results of a second GWAS were reported, in which the authors found eight SNPs involving seven genes to be associated with POP at a genome-wide significant level (p < 5.00 × 10 −8 ) [13]. As is the case in many other fields of medicine, the translational step from genetic data towards clinically useful information-i.e., diagnostic biomarkers and disease-modifying treatments-still needs to be made in the field of POP. Therefore, in this study, we conducted the first exome chip study of POP-which suggested multiple novel candidate genes-and integrated the exome chip results with other available genetic and expression data to generate and corroborate a molecular landscape of POP. This landscape provides detailed, novel insights into the biological processes that underpin pelvic floor function and dysfunction.

Exome Chip Study
We conducted an exome chip study with 526 POP cases and 960 control women (without POP). The characteristics of all participants in the exome chip study are presented in Table 1. Compared to cases, controls were on average 5 years older and had a higher ratio of postmenopausal females. Since the prevalence of POP increases with age and is probably associated with postmenopausal state, this means that there would likely have been more POP-associated SNVs identified if the case group would have had a similar average age and a similar percentage of postmenopausal women than the control group. Therefore, the results of our study may be underestimating the real genetic differences. Since the prevalence of POP increases with age and is probably associated with a postmenopausal state, the results of our study may be underestimating the real genetic differences. We identified 54 statistically significant-i.e., with Bonferroni-corrected p < 0.05, uncorrected p < 1.48 × 10 −6 -exonic and (potentially) deleterious single nucleotide variants (SNVs), including 52 missense mutations, 1 splice site mutation, and 1 stop-gain mutation. The list of these SNVs and the genes they affect is provided in Table 2. For all 54 SNVs, a combined annotation-dependent depletion (CADD) score was available-a single pathogenicity score that combines 63 distinct computational prediction methods [14], including the Grantham score for missense mutations [15]-and 18 SNVs had a CADD score > 20, while 12 had a CADD score > 10 but < 20 (Table 2). These SNVs are predicted to be among the 1 and 10% most deleterious of all possible substitutions in the human genome, respectively [14,16]. This implies that a considerable proportion of the SNVs that we identified have a damaging effect on protein function.  Abbreviations: NA, not applicable; POP, pelvic organ prolapse; SNV(s), single nucleotide variant(s). a A MAF of '0' indicates that none of the POP cases or controls were found to carry the minor allele of the SNV; b the Grantham score (GS) for missense mutation-type SNVs ranges from 0 to 215 and a higher score indicates a stronger negative effect of the amino acid change on protein structure and function [15]; c the combined annotation-dependent depletion (CADD) tool combines 63 distinct computational prediction methods-including the GS for missense mutations-and calculates a single pathogenicity score [14]. SNVs with CADD scores > 20 and > 10 indicate that these are predicted to be among the 1 and 10% most deleterious of all possible substitutions in the human genome, respectively [14,16]. Within our set of differentially expressed genes in an independent POP cohort (see Methods), the expression of six of the genes from the exome chip was different: the expression of ANKRD49 was downregulated (FC = −1.19; corrected p = 7.08 × 10 −3 ) d , the expression of ZNF700   which is initiated when epithelial cells start producing less CDH1 and more of the mesenchymal marker N-cadherin (CDH2) (not shown). This results in loss of epithelial cell-cell adhesion and the epithelial cells gradually transforming into mesenchymal cells, which then further differentiate into fibroblasts. EMT involves the key transcription factor SMAD3 and is negatively regulated by three sex hormone-bound receptors-the estrogen receptors ESR1 and ESR2 and the progesterone receptor (PGR)-that each upregulate CDH1 expression and hence prevent EMT. Conversely, TGFB1 induces EMT as it downregulates CDH1 expression.
(2) Proteins of the major histocompatibility complex II (MHC II)-i.e., HLA-DQA1 and HLA-DQB1-are expressed in epithelial cells of the female reproductive tract, where they are involved in regulating the activity of the immune response through presenting foreign antigens to circulating T lymphocytes (not shown). TGFB1 regulates the expression of these two HLA proteins, in part through upregulating the expression of beta-2microglobulin (B2M), an EMT-inducing extracellular protein that mediates antigen presentation.
(3) Further, the ECM provides support to epithelial cells and fibroblasts, and modulation of the ECM is essential for many pathophysiological processes such as tissue growth, wound healing, and fibrosis. The ECM is composed of different molecules, including elastin (ELN) and collagen (COL) fibers, as well as proteins that crosslink or regulate the degradation of ELN and COL (e.g., LOX and LOXL1), various matrix metalloproteinases (MMPs) and their regulators, TIMP1 and TIMP2. Other important ECM proteins are ECM1, FBLN5, LAIR2, and LAMC1. TGFB1 is involved in ECM remodeling through regulating the expression of most ECM proteins in the landscape. (4) Since fibroblasts are responsible for the synthesis and secretion of the main ECM components through SMAD3-which, as indicated above, also plays a major role in EMT-fibroblast survival and apoptosis (and hence proper functioning) is another important landscape process. TRAIL, a cytokine of the tumor necrosis factor (TNF) family, regulates fibroblast survival and apoptosis through its receptors and signaling involving pro-apoptotic proteins such as NFKB and CASP3. As with the three other main processes, TGFB1 is a key regulator, e.g., through regulating the expression of TRAIL-R1 and CASP3 and activating NFKB. In addition, extracellular AGE (or advanced glycation end product) molecules binding to their receptor RAGE results in NFKB activation, CASP3 inhibition, and translocation of SMAD3 to the nucleus, which results in the AGE-RAGE complex stimulating the apoptosis of fibroblasts. Lastly, nuclear proteins such as HIF1A and PARP1 are also important regulators of fibroblast survival and apoptosis.

Molecular Landscape of POP
A molecular landscape was built in which 26 of the 54 proteins encoded by the genes from the exome chip study (Table 2) interact functionally, together with 41 other POP candidates as well as 2 POP-implicated microRNAs (Supplementary Table S1). The molecular landscape of POP is shown in Figure 1 and includes 10 proteins/molecules (indicated in white) that have not been directly linked to POP (yet) but show important functional interactions within the landscape.
The molecular POP landscape is located in epithelial cells, fibroblasts, and the surrounding extracellular matrix (ECM) of the female urogenital tract. Four main biological processes operate within the landscape, i.e., epithelial-mesenchymal transition (EMT), immune response activation, modulation of the ECM, and fibroblast survival and apoptosis. These processes interact with each other and are regulated through signaling cascades involving female sex hormones and the cytokine transforming growth factor beta 1 (TGFB1). In the legend of Figure 1, a succinct description of these four processes is given, and in the Supplementary Information, a detailed description of the molecular signaling cascades and protein interactions in the landscape is provided.

Enrichment Analyses of Differential Expression Data from an Independent POP Patient Cohort
The upstream regulator analysis of the 618 genes that were differentially expressed between POP and non-POP tissues in the independent cohort of 12 premenopausal women with POP (see Materials and Methods) at a corrected p < 0,01 (Supplementary Tables S2 and S3) yielded five regulators that are also present in the POP landscape: TGFB1, SMAD3, PGR and ESR2 are predicted to be activated while COL18A1 is predicted to be inhibited (Table 3).
Moreover, the targets of TGFB1 are most significantly enriched within the set of 618 genes. We identified five functional gene categories that are present in the POP landscape and are enriched within the 618 differentially expressed genes: three categories relating to (epithelial) cell death-i.e., 'apoptosis', 'necrosis of epithelial tissue' and 'cell death of epithelial cells'-are predicted to be inhibited, while two categories relating to (fibroblast) survival-i.e., 'cell survival' and 'differentiation of fibroblasts'-are predicted to be activated (Table 4). We also identified five genes that encode proteins within the POP landscape and are differentially expressed between POP and non-POP tissues in the independent POP patient cohort, i.e., AKR1C1, PLAUR, and VIM are upregulated while FBLN5 and RAD52 are downregulated (Supplementary Table S4). Table 3. Upstream regulator analysis of the 618 mRNAs that were differentially expressed (corrected p < 0.01) between POP and non-POP tissues in 12 POP patients (see Materials and Methods). The 5 upstream regulators that are present in the POP landscape are listed with a z-score ≥ 2.00 or ≤ −2.00, indicating that they are predicted to be activated or inhibited, respectively. For each regulator, the p-value-referring to how enriched the targets of each upstream regulator are in the differentially expressed mRNAs-and the number of target genes are also listed.

Metformin Downregulates the Expression of TGFB1 Target Genes in Vaginal Fibroblasts from POP Patients
As TGFB1 is the 'master regulator' of the POP landscape and its activation results in many landscape processes that are involved in POP, we tried to confirm metformin-an oral anti-diabetic drug that counteracts TGFB1-as a drug target 'lead' for POP by testing its effect on the expression of selected genes from the landscape. Since metformin suppresses TGFB1 responses specifically in fibroblasts cell lines [17], we used human vaginal fibroblasts from four independent POP patients to see if the drug could have the same beneficial effect in POP. As shown in Figure 2, we found that metformin (negatively) affects the TGFB1induced expression of three ECM proteins in these fibroblasts: COL3A1 (p = 0.0136) and ELN (p = 0.039) were significantly downregulated, while there was a tendency towards lower expression of COL1A1 (p = 0.0961).

Discussion
We performed the first exome chip study of POP and integrated the results with the available genetic and expression data into a protein interaction landscape of the disease. This molecular landscape contains the proteins encoded by 26 of the 54 genes affected by SNVs from the exome chip study, as well as 41 other POP candidate molecules and two POP-implicated microRNAs. The POP landscape reveals four main biological processes that operate in epithelial cells, fibroblasts, and the surrounding ECM of the female urogenital tract: EMT, immune response activation, ECM modulation, and fibroblast survival and apoptosis. During pregnancy and (vaginal) delivery, similar changes in immune response activation and ECM modulation than in the POP landscape occur [20,21], which is

Discussion
We performed the first exome chip study of POP and integrated the results with the available genetic and expression data into a protein interaction landscape of the disease. This molecular landscape contains the proteins encoded by 26 of the 54 genes affected by SNVs from the exome chip study, as well as 41 other POP candidate molecules and two POP-implicated microRNAs. The POP landscape reveals four main biological processes that operate in epithelial cells, fibroblasts, and the surrounding ECM of the female urogenital tract: EMT, immune response activation, ECM modulation, and fibroblast survival and apoptosis. During pregnancy and (vaginal) delivery, similar changes in immune response activation and ECM modulation than in the POP landscape occur [20,21], which is interest-ing as increased parity and vaginal delivery are established environmental risk factors for developing POP [5].
Of note, none of the "classical" POP candidate genes-identified through candidate gene association studies and mainly encoding ECM proteins-or the genes from two previous GWASs of POP (including the landscape genes COL18A1, FBLN3, NOP56, TBX5, and WNT4) contained significant SNVs in our study.
Our landscape does not assume or represent a "sequence of events", i.e., a number of biological processes and signaling cascades that occur in a spatially and temporally distinct order. Instead, a deficit in any or a combination of the four main landscape processes or their constituting signaling cascades, possibly aggravated through environmental risk factors and caused by variants in one or more genes, can lead to POP in various degrees of severity. This is reflected by the clinical reality of POP as a disease with a multifactorial etiology and variable treatment response.
Multiple signaling cascades in the landscape are regulated by female sex hormones (i.e., estrogen in its active form, estradiol, and progesterone). These hormones indirectly affect oxidative stress levels in epithelial cells and fibroblasts of the urogenital tract, as progesterone-activated PGR [22] and estradiol-activated ESR2 [23] upregulate the expression of the enzyme GGCT. Together with three other enzymes from the landscape-GPX1, GSTP1, and METAP1-GGCT is involved in the metabolism of glutathione [24], an important antioxidant that prevents oxidative stress-induced cellular damage by reactive oxygen species (ROS) [25]. Interestingly, ROS production by fibroblasts of women without POP was found to be increased through mechanical stress [26], while oxidative stress biomarkers were also increased in the uterosacral ligaments of POP patients [27]. ROS can also directly induce EMT [28], corroborating the role of oxidative stress in POP. Although sex hormones influence many landscape cascades, their role in POP etiology and, hence, their clinical relevance, is not completely understood. Increasing age and postmenopausal state are risk factors for POP [4,5], and the hypoestrogenic state after the menopause could be a factor in disease onset. Both lower serum levels of estradiol [29][30][31] and lower serum and urine levels of progesterone (metabolites) [29,32] have been associated with an increased risk of POP. Therefore, hormone replacement therapy (HRT) in postmenopausal women could have a favorable effect on disease progression or even prevent the disease from developing. However, studies on the effect of HRT in POP patients show conflicting results, as negative [33], protective [34], or no effects [35][36][37] have all been reported. Given the above, we suggest that future studies on the effect of HRT on POP take into consideration confounding factors such as endogenous sex hormone levels, menopause duration, or known associated genetic variations.
In addition to sex hormones and oxidative stress, the cytokine TGFB1 is an important landscape protein, as it represents the 'master regulator' that, when activated, affects and modulates all four biological processes in the landscape. Notably, (activated) TGFB1 is involved in ECM remodeling and regulating fibroblast function. Moreover, TGFB1 itself has an important role in POP etiology, and TGFB1 is both differentially expressed [31,38] and could have a functional effect-e.g., through affecting fibroblast proliferation [39]-in POP patients.
We have been able to corroborate our findings by analyzing differential gene expression data from an independent POP cohort. Our analysis of these data revealed that four 'upstream regulators'-the most significant regulator TGFB1 and three other proteins that are important in EMT, i.e., SMAD3 and the sex hormone receptors PGR and ESR2-are predicted to be activated in POP tissue, while COL18A1, a collagen with multiple roles in ECM modulation, is predicted to be inhibited. These results are in line with the landscape, in which EMT and ECM modulation are two of the four main processes. Moreover, five genes encoding proteins in the landscape-including AKR1C1 and RAD52, two genes from the exome chip study-are differentially expressed in the independent cohort, providing a direct validation of these genes as POP candidates. Further, our findings from the enrich-ment analysis of functional gene categories indicate that both epithelial cells and fibroblasts are affected in POP.
As the 'master regulator' of the landscape, TGFB1 would be an interesting target for developing novel treatments for POP. More specifically, drugs that lower TGFB1 activity could have beneficial effects and eventually be used in a clinical setting. However, as a result of its diverse effects on multiple signaling cascades in the body, drugs (negatively) targeting the expression or function of the genes encoding TGFB1 or its receptors are likely to be associated with multiple side effects. Hence, it would be better to test drugs that selectively target individual signaling pathways downstream of TGFB1. Such a novel treatment for POP may be metformin, an oral anti-diabetic drug that affects multiple genes/proteins, and processes in the landscape. First, metformin reduces the expression of TGFB1 [40,41] and decreases the TGFB1-promoted loss of CDH1 [42,43], which results in less EMT. Metformin also regulates female sex hormone signaling as it suppresses the expression of miR-221 [44,45] and miR-222 [44], two POP-implicated microRNAs that downregulate ESR1 expression. Metformin inhibits the TGFB1-induced synthesis of type 1 collagen by fibroblasts through suppressing SMAD3 phosphorylation [46][47][48] (a protein that has been found to be upregulated in POP tissues [49]). Furthermore, metformin may have a (more) directly beneficial effect on cells from the urogenital tract (and hence on POP) as it was found to reduce estradiol-induced EMT [50], inflammation [51], and age-associated dysfunction [52] of uterine cells. Therefore, since metformin also specifically suppresses TGFB1 responses in fibroblast cell lines [17], we used fibroblasts from patients to assess if the drug could affect TGFB1 signaling in a direction that might have a beneficial effect on POP (or POP development). We indeed showed that metformin counteracts the molecular mechanisms downstream of activated TGFB1 that negatively affect POP, through downregulating the expression of the key ECM (and landscape) proteins COL1A1, COL3A1, and ELN. These in vitro findings add evidential weight to metformin as a novel POP treatment, and further (pre)clinical studies are needed to fully develop it as a repurposed drug.
The current results should be viewed in light of some strengths and limitations. The main limitations of our study are that we have no independent replication of the exome chip results and we have only initial in vitro findings for the beneficial effect of one drug that could be repurposed for POP treatment. However, our findings are strengthened by the enrichment analyses of differential expression data from an independent POP cohort. In addition, our findings for metformin by no means show clinical effectiveness but could serve as the starting point for further studies on this and other drugs that could be repurposed.

Participant Characteristics and Quality Control
A case-control design was used for the exome chip study. Data on 242,901 genetic markers and 2478 individuals were available prior to the implementation of Quality Control (QC), which was implemented using PLINK software, version 1.9. Of the 2478 individuals, 548 were classified as cases, while the remaining 1930 individuals were considered to be controls. Given the nature of the studied phenotype, all cases were female, while the controls consisted of 977 females and 953 males. Prior to QC, the total genotype calling rate was 0.999053. A total of 12 individuals (11 cases and 1 control) were removed for a low per-individual calling rate (<0.95), while the same per-marker calling rate criterion (<0.95) led to the removal of 82 markers. Further QC steps in the combined case/control population included checks of markers for low minor allele frequency (MAF) (MAF < 0.005, failed by 199,851 markers), deviation from Hardy-Weinberg equilibrium (HWE) (p < 10 −4 , 280 markers), differential missingness with respect to phenotype (via--test-missing, p < 10 −3 , 88 markers), and haplotype (via--test-mishap, p < 10 −6 , 14 markers). Furthermore, individuals were checked for sex discordance (via--check-sex, 10 individuals), individual duplicity or high-relatedness (via--genome, PI_HAT ≥ 0.9 [53], 7 individuals), and heterozygosity (via--het, with conf. int. based on MAD, 13 individuals). Heterozygous haploid genotypes and non-missing non-male Y chromosome genotypes were set to missing values. Lastly, all markers that are not located within a gene were removed (8732 markers). After all the QC steps and removing all the male controls-to allow for a 'clean' comparison with the allfemale POP cases-there were 33,854 markers and 1486 females (526 cases and 960 controls) left for the exome chip analysis. The case group of 526 women presented with symptoms and signs of POP (including women with recurrent POP and previous POP surgery) at the Department of Obstetrics and Gynecology of the Radboud University Medical Center (Radboudumc), Nijmegen, The Netherlands, in January 2007-February 2011. Signs of POP were defined as stage II or more, according to the Pelvic Organ Prolapse Quantification (POP-Q) system [54], i.e., the descent of the leading edge of the prolapse 1 cm above the hymenal remnants or beyond. POP had already been surgically treated in some women at the time of inclusion. Exclusion criteria were a poor understanding of the Dutch language and/or genetic diseases with an increased risk of POP (e.g., Ehlers-Danlos and Marfan syndromes, and myotonic dystrophy type 1). For POP cases, we collected data on previous surgery, POP-Q stage, as well as POP complaints measured with the validated Dutch translation of the Urogenital Distress Inventory (UDI) [55]. If more than one measurement was available over time, the most severe measurement was used. The UDI records the presence or absence of specific symptoms, and associated bother on a four-point Likert scale. The UDI scores are transformed to a continuous scale, ranging from 0 (symptom absent or no bother) to 100 (a very bothersome symptom). According to this scale, women without symptoms and women with symptoms but without any bother are scored equally at 0. The exome chip study was approved by the local institutional review board in an amendment of CMO numbers 2007/043 and 2010/071. As the control group, we used the exome chip data available from 960 women from the Nijmegen Biomedical Study (NBS), a general population-based survey conducted by Radboudumc's Departments for Health Evidence and Department of Laboratory Medicine [56]. For the cases and controls, data on age, parity, and BMI were collected.

Laboratory Analysis and Genotyping
Blood samples were drawn from all cases and controls. Genomic DNA was isolated from peripheral blood using the Nucleospin Blood Kit (Macherey-Nagel ® , Düren, Germany). Both case and control samples were genotyped (Illumina Exome BeadChip v1.1, Illumina ® San Diego, CA, USA) and called at the Human Genome Facility and the department of Epidemiology, Erasmus MC, Rotterdam. Calling was performed in GenomeStudio using the default cluster files provided by Illumina (HumanExome-12v1-1_A.egt), followed by the calling of non-called variants using zCall [57].

Population Structure
Since principal component analysis uncovered evidence of population structure in the data, we performed clustering analysis using the R function hclust() with method ward.D [58]. The dissimilarity structure was generated by dist() using the first 5 components produced by PLINK's multidimensional scaling (MDS) analysis (via--mds-plot 10) to estimate continuous axes of ancestry. Using the output of the R function NbClust(), the desired number of clusters was set to 3. The case/control counts in the resulting three clusters were 96/182,124/195 and 306/583.

Statistical Analysis
Using the quality controlled clustered data, we analyzed the 33,854 markers for association with POP using the Cochran-Mantel-Haenszel (CMH) test's asymptotic p-values via PLINK's --mh option at the Bonferroni-corrected level of significance (0.05/33,854 = 1.48 × 10 −6 ).

Molecular Landscape Building
Applying the approach that we previously used for multiple neurodevelopmental and neurological disorders [59][60][61], we then built a molecular landscape of POP. First, we searched the literature for the (putative) function of all candidate genes implicated through the exome chip study as well as other POP candidates implicated via other evidence, including (genome-wide) genetic association studies, miRNA/mRNA/protein expression studies, and/or genetic animal model studies (Supplementary Table S1). For this literature analysis, we used the UniProt Protein Knowledge Base (http://www.uniprot.org/uniprot accessed on 21 June 2022) to gather basic information on the function of all candidate genes and their encoded proteins. Based on the biological processes that have already been suggested to be impaired in POP, we then searched PubMed (http://www.ncbi.nlm.nih. gov/sites/entrez accessed on 21 June 2022) for all genes/proteins and using the search terms "collagen", "connective tissue", "fibrosis", "pelvic" and "pelvic organ". Further, we searched PubMed for all functional interactions among all candidate proteins. Lastly, based on all the gathered information, we generated a molecular landscape of protein interactions. The figure of the landscape was made using the drawing program Serif DrawPlus version 4.0 (www.serif.com accessed on 21 June 2022).

Enrichment Analyses of Differential Gene Expression Data from an Independent POP Patient Cohort
In order to corroborate our findings, we conducted enrichment analyses of genes that were differentially expressed between POP and non-POP tissues in POP patients. For this, we used microarray data that had been previously deposited by authors MHK and AMRZ on the GEO website (http://www.ncbi.nlm.nih.gov/gov/geo/ accessed on 21 September 2022). The data set GSE53868 contains microarray-derived gene expression data from 12 premenopausal women with POP, and the data are available for POP (prolapsed tissue) and non-POP (peri-cervical region) tissues from the anterior vaginal wall of each patient. To identify expression changes between POP and non-POP tissues, we performed differential gene expression analysis using R statistical software version 4.0 and the R/Bioconductor [62] package limma [63]. Since there were two arrays per individual patient, we used the "paired samples" design, allowing for correlation analyses within patients. The Benjamini-Hochberg (BH) method was used to correct for multiple testing, and only protein-coding genes with an adjusted p-value <0.01, independent of magnitude of change, were considered as differentially expressed in pair-wise comparisons between POP and non-POP sites, and were used in the subsequent enrichment analyses.
First, we used the Ingenuity Pathway Analysis (IPA) software (winter release 2022) (Qiagen Bioinformatics, https://www.qiagenbioinformatics.com/products/ingenuity-pathwayanalysis/ accessed on 21 December 2022) to perform gene enrichment analyses of the 618 unique genes that were differentially expressed between the POP and non-POP tissues of 12 POP patients at a BH-corrected p-value < 0.01 (Supplementary Table S2). We then conducted an upstream regulator analysis to identify the proteins that are both present in the POP landscape and regulate the expression of multiple targets within the set of 618 genes. Based on the "Ingenuity Knowledge Base", a repository of data from publicly accessible databases and data that are manually curated by systematically reviewing published literature, IPA can generate a list of 'upstream regulators', i.e., proteins or compounds that regulate the expression of multiple target genes, from an input list. The program also calculates a z-score based on the expression changes of the input genes, which is a measure of the directionality of the upstream regulator. In this respect, a z-score ≥2.00 and ≤−2.00 (reflecting activation and inhibition of the upstream regulator-dependent effects on target gene expression, respectively) is considered to be significant. IPA also calculates a p-Value that reflects how enriched the targets of each upstream regulator are in the differentially expressed genes. In addition to generating a list of upstream regulators, IPA assigns genes and their corresponding proteins to functional categories and, as with the upstream regulators, we identified those functional categories that are both present in the POP landscape and are enriched in the set of 618 genes. IPA also generated z-scores-reflecting whether the functional gene category was activated or inhibited-as well as BH-corrected P-values of enrichment for each category. Further, for all genes in the POP landscape, we searched the list of all differentially expressed genes (Supplementary Table S3) to ascertain whether their expression was (significantly) upregulated or downregulated (Supplementary Table S4).

In Vitro Experiments and Statistical Analysis
In this study, we used primary human vaginal fibroblasts (passages 6-7) from the prolapsed anterior vaginal wall of four premenopausal women with cystocele that were previously isolated from our group [64]. Full thickness biopsies (>1 cm 2 ) from the anterior vaginal wall were retrieved from each of the women. Biopsies were collected in PBS at 4 • C, and within 24 h, primary human vaginal fibroblasts were isolated, as previously described [65]. Experiments were performed in an incubator at 37 • C, 95% humidity, and 5% CO 2 . A cell density of 10,000 cells/cm 2 were seeded on 6-well plates and cultured for 24 h with the cell culture medium DMEM (Gibco-Life technologies, Paisley, UK) supplemented with 10% fetal bovine serum (FBS; Sigma-Aldrich, St. Louis, MO, USA), 100 µg/mL streptomycin, 100 U/mL penicillin, and 250 µg/mL amphotericin-B (Sigma). Subsequently, the media was refreshed with cell culture media supplemented with 0.1 ng/mL of human TGF beta 1 recombinant protein (eBioscience, Thermo Fisher Scientific, Waltham, CA, USA). After 24 h of culture, cells were stimulated (or not) with 2 mM of metformin hydrochloride (Tocris Bioscience, Bristol, UK) for an additional 24 h. At the end of the experiment, total RNA was isolated using TRIzol and following the manufacturer's instructions (Invitrogen). Concentration and purity of the RNA were determined with a Nanodrop-1000 spectrophotometer (Thermo Scientific, Waltham, CA, USA). Subsequently, 2 ug DNase-I-treated total RNA was used to make cDNA using random hexamer primers (Roche) and SuperScript II Reverse Transcriptase (Invitrogen). Real-time PCR (RT-PCR) was performed with a LightCycler 480 SYBR Green I master mix and in a LightCycler LC480 device (Roche) with the following amplification conditions: 5 min at 95 • C, followed by 45 cycles of 10 sec at 95 • C, 20 sec at 60 • C, and 20 sec at 72 • C. Crossing-point (Cp) values were determined using the LightCycler480 SW, version 1.5 software (Roche). The gene expression levels for three selected genes from the landscape-the collagen genes COL1A1 and COL3A1, and elastin (ELN)-were normalized to the housekeeping genes YWHAZ and HPRT1 and were calculated for the vaginal fibroblasts from each of the women with POP separately and according to the mathematical model for relative quantification in RT-PCR described by Pfaffl [19]. COL1A1, COL3A1 and ELN were selected as they encode key, TGFB1-induced components of the ECM surrounding fibroblasts in the anterior vaginal wall that regulate the mechanical strength of this wall, and as such, are dysregulated in POP. Further, to find the right housekeeping genes, we first performed a pilot experiment in which we tested which common housekeeping genes show stable expression in vaginal fibroblasts regardless of the conditions these cells are in [66]. YWHAZ and HPRT1 turned out to be the most stable and were therefore chosen as reference housekeeping genes. The primers for the three landscape genes and the two housekeeping genes were acquired from Life Technologies and are listed in Supplementary Table S5. The results were calculated as the percentage of downregulated gene expression by metformin in TGFB1-stimulated fibroblasts and are presented in a scatter plot. A one-sample t-test compared to 100 was used to test the effect of metformin on TGFβ-stimulated cells with the software Prism version 5.03 (GraphPad Software Inc., San Diego, CA, USA). We also applied Grubbs' test [18] for outlier correction.

Conclusions
Based on the results of our exome chip study and available literature data, we built a molecular landscape of POP that provides novel insights into the biological processes underlying the disease and that we could corroborate through analyzing gene expression data from an independent cohort. In addition, the POP landscape provides interesting leads for existing drugs such as metformin that could be further developed as repurposed compounds in selected patient groups.
Supplementary Materials: The supporting information can be downloaded at: https://www.mdpi. com/article/10.3390/ijms24076087/s1. Funding: The exome chip data for the patients was funded by the Department of Obstetrics and Gynecology of the Radboudumc (Nijmegen, The Netherlands), and NBS controls were generated in a research project that was financially supported by BBMRI-NL, a research infrastructure financed by the Dutch government (NWO 184.021.007) (https://www.nwo.nl/ accessed on 21 December 2022). Co-author KAB was supported by NIH grant RO1 HD061821 from the Eunice Kennedy Shriver National Institute of Child Health and Human Development (https://www.nih.gov/ accessed on 21 December 2022). Co-authors KBK, AMRZ, JW, WDW, WMP, EO, and GP were supported by the European Regional Development Fund (ERDF) (https://ec.europa.eu/regional_policy/en/funding/ erdf/ accessed on 21 December 2022) under grant agreement no. PROJ00787 (DIABIP). The funders had no role in study design, data collection and analysis, the decision to publish, or the preparation of the manuscript.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki. The exome chip study was approved by the local institutional review board in an amendment of CMO numbers 2007/043 and 2010/071. As the control group, we used the exome chip data available from 960 women from the Nijmegen Biomedical Study (NBS, a general populationbased survey conducted by Radboudumc's Departments for Health Evidence and Department of Laboratory Medicine [56]).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
All the key data that support the findings of this study are available in the main text or the Supplementary Materials. Other, more extensive data are available from the corresponding author upon reasonable request.