Disrupted PGR-B and ESR1 signaling underlies defective decidualization linked to severe preeclampsia

  1. Tamara Garrido-Gomez  Is a corresponding author
  2. Nerea Castillo-Marco
  3. Mónica Clemente-Ciscar
  4. Teresa Cordero
  5. Irene Muñoz-Blat
  6. Alicia Amadoz
  7. Jorge Jimenez-Almazan
  8. Rogelio Monfort-Ortiz
  9. Reyes Climent
  10. Alfredo Perales-Marin
  11. Carlos Simon  Is a corresponding author
  1. Igenomix Foundation, INCLIVA, Spain
  2. Igenomix, Spain
  3. Department of Obstetrics and Gynecology, University and Polytechnic La Fe Hospital, Spain
  4. Department of Obstetrics and Gynecology, School of Medicine, Valencia University, Spain
  5. Obstetrics & Gynecology, BIDMC Harvard University, United States

Abstract

Background:

Decidualization of the uterine mucosa drives the maternal adaptation to invasion by the placenta. Appropriate depth of placental invasion is needed to support a healthy pregnancy; shallow invasion is associated with the development of severe preeclampsia (sPE). Maternal contribution to sPE through failed decidualization is an important determinant of placental phenotype. However, the molecular mechanism underlying the in vivo defect linking decidualization to sPE is unknown.

Methods:

Global RNA sequencing was applied to obtain the transcriptomic profile of endometrial biopsies collected from nonpregnant women who suffer sPE in a previous pregnancy and women who did not develop this condition. Samples were randomized in two cohorts, the training and the test set, to identify the fingerprinting encoding defective decidualization in sPE and its subsequent validation. Gene Ontology enrichment and an interaction network were performed to deepen in pathways impaired by genetic dysregulation in sPE. Finally, the main modulators of decidualization, estrogen receptor 1 (ESR1) and progesterone receptor B (PGR-B), were assessed at the level of gene expression and protein abundance.

Results:

Here, we discover the footprint encoding this decidualization defect comprising 120 genes—using global gene expression profiling in decidua from women who developed sPE in a previous pregnancy. This signature allowed us to effectively segregate samples into sPE and control groups. ESR1 and PGR were highly interconnected with the dynamic network of the defective decidualization fingerprint. ESR1 and PGR-B gene expression and protein abundance were remarkably disrupted in sPE.

Conclusions:

Thus, the transcriptomic signature of impaired decidualization implicates dysregulated hormonal signaling in the decidual endometria in women who developed sPE. These findings reveal a potential footprint that could be leveraged for a preconception or early prenatal screening of sPE risk, thus improving prevention and early treatments.

Funding:

This work has been supported by the grant PI19/01659 (MCIU/AEI/FEDER, UE) from the Spanish Carlos III Institute awarded to TGG. NCM was supported by the PhD program FDGENT/2019/008 from the Spanish Generalitat Valenciana. IMB was supported by the PhD program PRE2019-090770 and funding was provided by the grant RTI2018-094946-B-100 (MCIU/AEI/FEDER, UE) from the Spanish Ministry of Science and Innovation with CS as principal investigator. This research was funded partially by Igenomix S.L.

Introduction

Preeclampsia (PE) is a severe complication of late pregnancy and is the second leading cause of maternal mortality in the US, affecting 8% of first-time pregnancies (Gifford, 2000). PE is characterized by the onset of hypertension, proteinuria, and other signs of maternal vascular damage that contributes to maternal and neonatal mortality and morbidity (Gifford, 2000). Severe preeclampsia (sPE) is diagnosed based on elevated blood pressure (systolic ≥160 or diastolic of ≥100 mm Hg) or thrombocytopenia, impaired liver function, progressive renal insufficiency, pulmonary edema, or the onset of cerebral or visual disturbances (Gynecologists ACoOa, Pregnancy TFoHi, 2013). sPE is a placental insufficiency syndrome mediated by early-deficient extravillous trophoblast (EVT) invasion of uterine decidua and spiral arterioles, leading to incomplete endovascular invasion and altered uteroplacental perfusion (Roberts and Cooper, 2001; Brosens et al., 2019; Staff et al., 2020). Why shallow EVTs invasion occurs, however, remains to be determined (Fisher, 2015).

Pregnancy health is dictated by the embryo, placenta, and the quality of the maternal decidua, where EVTs invasion and remodeling of maternal spiral arteries occur (Norwitz et al., 2001; Cha et al., 2012). Accumulated evidence suggests that the contribution of the decidua to the etiology of PE (Rabaglino et al., 2015), sPE (Garrido-Gomez et al., 2017; Garrido-Gomez et al., 2020; Garrido-Gómez et al., 2020), and placenta accreta (McNally et al., 2020) is significant, and cellular signaling in the decidua may determine whether these conditions develop. Decidualization is the remodeling of the endometrium initiated after ovulation that is necessary for adequate trophoblast invasion and subsequent placentation (Gellersen and Brosens, 2014). Defective decidualization (DD) entails the inability of the endometrial compartment to undertake tissue differentiation, leading to aberrations in placentation and compromising pregnancy health (Garrido-Gómez et al., 2020).

In humans and other great apes, the formation of the decidua is a conceptus-independent process driven by progesterone and the second messenger cyclic adenosine monophosphate (Brar et al., 1997) that stimulates synthesis of a complex network of intracellular and secreted proteins through progesterone receptor (PR) activation. Endometrial decidualization involves secretory transformation of uterine glands (Kelleher et al., 2018), influx of specialized immune cells, vascular remodeling, and morphological (Dunn et al., 2003; Ramathal et al., 2010), biochemical (Giudice et al., 1998; Jabbour and Critchley, 2001), and transcriptional reprogramming of the endometrial stromal compartment (Wang et al., 2020). We recently characterized the transcriptomics of human decidualization at single-cell resolution from secretory endometrial samples and showed that the process is initiated gradually after ovulation, with a direct interplay between stromal fibroblasts and lymphocytes (Wang et al., 2020). However, most knowledge on decidual function in health and disease comes from in vitro model systems (Ng et al., 2020; Zhou et al., 2013; Ganeff et al., 2009).

In the present study, we aimed to discern the preconception decidual transcriptomic signature associated with in vivo DD. We performed a comparative global transcriptional profiling of endometrium in women who developed sPE in a previous pregnancy. Initially, we identified 593 genes differentially expressed in sPE compared to control cases. Then, molecular DD fingerprint of 120 genes associated with the development of sPE was defined and evaluated as a diagnostic tool in an independent cohort of samples. Finally, we identified PR and estrogen receptor 1 (ER1) as targets of the DD fingerprint genes. Our findings indicate that an endometrial transcriptomic signature persists years after the affected pregnancy. This signature may be leveraged for a preconception or early prenatal screening strategy in assessing sPE risk and may inform the development of sPE therapies.

Methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Biological sample (Homo sapiens)Endometrial biopsiesUniversity and Polytechnic La Fe Hospital (Valencia, Spain)Freshly isolated from human donors
AntibodyAnti-progesterone receptor antibody [YR85] (rabbit monoclonal anti-human)AbcamCat: AB32085RRID:AB_777452Dilution: (1:50)
AntibodyAnti-estrogen receptor alpha antibody (mouse monoclonal antibody)Santa CruzCat: sc-8002RRID:AB_627558Dilution: (1:50)
AntibodyGoat anti-rabbit IgG H&L (Alexa Fluor 488) (goat polyclonal)AbcamCat: ab150077RRID:AB_2630356Dilution: (1:1000)
AntibodyGoat anti-mouse IgG (H + L) Cro Alexa Fluor 488 (goat polyclonal)InvitrogenCat: A-11001RRID:AB_2534069Dilution: (1:1000)
Sequence-based reagentRT-qPCR primersThis paperSupplementary file 3
Commercial assay or kitQIAsymphony RNA KitQiagen931636Global RNA-seq library preparation
Commercial assay or kitIllumina TruSeq Stranded mRNA sample prep kitIllumina20020595Global RNA-seq library preparation
Commercial assay or kitKapa SYBR fast qPCR kitKapa Biosystems IncKK4602Global RNA-seq library preparation
Commercial assay or kitTruSeq RNA CD Index Plate (96 indexes, 96 samples)Illumina20019792RNA sequencing
Commercial assay or kitNextSeq 500/550 cartridge of 150 cyclesIlluminaFC-404-2002RNA sequencing
Commercial assay or kitSuperScript VILO cDNA Synthesis KitThermo Fisher Scientific11754250RT-qPCR. cDNA preparation
Software, algorithmSTARDobin et al., 2013URL: http://code.google.com/p/rna-star/RRID:SCR_004463RNA-seq analysisRead alignerVersion 2.4.2a
Software, algorithmFastQCURL: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/RRID:SCR_014583RNA-seq analysisQuality of FASTQ file determinationVersion 0.11.2
Software, algorithmSAMtoolsLi et al., 2009URL: http://htslib.org/RRID:SCR_002105RNA-seq analysisSAM and BAM manipulation files Version 1.1
Software, algorithmHTSeqAnders et al., 2015URL: http://htseq.readthedocs.io/en/release_0.9.1/RRID:SCR_005514RNA-seq analysisTo count the number of reads per geneVersion 0.6.1p1
Software, algorithmBEDtoolsQuinlan and Hall, 2010URL: https://github.com/arq5x/bedtools2RRID:SCR_006646RNA-seq analysisTo obtain gene coverageVersion 2.17.0
Software, algorithmedgeRRobinson et al., 2010URL: http://bioconductor.org/packages/edgeR/RRID:SCR_012802RNA-seq analysisTo analyze differentially expressed genesVersion 3.24.3
Software, algorithmStringJensen et al., 2009URL: http://string.embl.de/RRID:SCR_005223Interaction Network.
Software, algorithmCytoscapeShannon et al., 2003URL: http://cytoscape.orgSCR_003032Interaction Network
Software, algorithmCytoHubbaChin et al., 2014URL: http://apps.cytoscape.org/apps/cytohubbaRRID:SCR_017677Interaction Network
OtherCustom scriptsURL: https://github.com/mclemente-igenomix/garrido_et_al_2021The specific script to run RNA-seq analysis

Study design

A total of 40 non-pregnant women who experienced a previous pregnancy were enrolled in this study for endometrial RNA-sequencing analysis. Endometrial samples were obtained for research purposes during late secretory phase in 24 women who had developed sPE in a previous pregnancy and in 16 women with no history of sPE with full term (n = 8) and preterm pregnancies (n = 8) as controls. sPE was clinically defined based on elevated blood pressure (systolic ≥160 or diastolic of ≥100 mm Hg) or thrombocytopenia, impaired liver function, progressive renal insufficiency, pulmonary edema, or the onset of cerebral or visual disturbances. Endometrial biopsies were processed to obtain RNA and then converted to cDNA for library generation to perform next-generation sequencing. The experimental design was based on a stratified random sampling with a 70:30 proportion in two cohorts: a training (n = 29) and validation (n = 11) set of samples. The training set of samples was analyzed by RNA-seq to identify the global transcriptomic profiling changes between control (n = 12) and sPE (n = 17) samples. Selection criteria were applied to define a transcriptomic fingerprinting associated with DD detected in sPE. Finally, targeted analysis of the DD signature was validated in the test set composed of controls (n = 4) and sPE (n = 7).

Human donors

Endometrial samples were collected from women aged 18–42 without any medical condition who had been pregnant 1–8 years earlier. All participants had regular menstrual cycles (26–32 days) with no underlying gynecological pathological conditions and had not received hormonal therapy in the 3 months preceding sample collection. After the inclusion criteria were applied, endometrial biopsies were obtained by pipelle catheter (Genetics Hamont-Achel, Belgium) under sterile conditions in the late secretory phase (cycle days 22–32). Specimens were kept in preservation solution until processing. Maternal and neonatal characteristics of women with sPE and controls are summarized in Supplementary file 1. Biological and technical variables for each donor were considered to discard confounding effects on the transcriptomic profile (Supplementary file 2). This study was approved by the Clinical Research Ethics Committee of University and Polytechnic La Fe Hospital (Valencia, Spain; 2011/0383), and written informed consent was obtained from all participants before tissue collection and all samples were anonymized.

RNA extraction

Total RNA from endometrial biopsies was isolated using QIAsymphony RNA kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol. RNA concentrations were quantified using a Multiskan GO spectrophotometer (Thermo Fisher Scientific, Waltham, USA) at a wavelength of 260 nm. Integrity of the total RNA samples was evaluated by the RNA integrity number (RIN) and DV200 metrics using an Agilent high-sensitivity RNA ScreenTape in a 4200 TapeStation system (Agilent Technologies Inc, Santa Clara, CA). Samples used for the global RNA-seq showed RIN values ranging from 4.9 to 9.2.

Global RNA-seq library preparation and transcriptome sequencing

cDNA libraries from total RNA samples (n = 40) were prepared using an Illumina TruSeq Stranded mRNA sample prep kit (Illumina, San Diego, CA) following a balanced batch-group design. 3 µg of total RNA were used as the RNA input according to the manufacturer’s protocol. mRNAs were isolated from the total RNAs by purifying the poly-A containing molecules using poly-T oligo attached to magnetic beads. The RNA fragmentation, first- and second-strand cDNA syntheses, end repair, single ‘A’ base addition, adaptor ligation, and PCR amplification were performed according to the manufacturer’s protocol. The average size of the cDNA libraries was approximately 350 bp (including the adapters). cDNA libraries were quantified using an Agilent D1000 ScreenTape in a 4200 TapeStation system (Agilent Technologies Inc). Libraries were normalized to 10 nM and pooled in equal volumes. The pool concentration was quantified by qPCR using the KAPA Library Quantification Kit (Kapa Biosystems Inc) before sequencing in a NextSeq 500/550 cartridge of 150 cycles (Illumina). Indexed and pooled samples were sequenced 150 bp paired-end reads by on the Illumina NextSeq 500/550 platform according to the Illumina protocol.

RNA-seq analysis

Reads were mapped to the hg19 human genome transcriptome using the STAR (version 2.4.2 a) read aligner (Dobin et al., 2013). FastQC (version 0.11.2) was used to determine the quality of FASTQ files. The manipulation of SAM and BAM files was done with the software SAMtools (version 1.1) (Li et al., 2009). To count the number of reads that could be assigned to each gene, we used HTSeq (version 0.6.1p1; Anders et al., 2015) and BEDtools software (version 2.17.0; Quinlan and Hall, 2010) to obtain gene coverage and work with bedFiles. Quality control filters in each program were used following the software package recommendations, and reads were filtered by mapping quality greater than 90%. Transcriptomic data were deposited in the Gene Expression Omnibus database (accession number GSE172381). The Bioconductor package edgeR (version 3.24.3; Robinson et al., 2010) was used to analyze differentially expressed genes (DEGs). The trimmed mean of M-values normalization method was applied to our gene expression values. The glmTreat function was used to find DEGs between groups. The p-value adjustment method was false discovery rate (FDR) with a cutoff of 0.05 (FDR < 0.05) and the fold-change (FC) threshold was 1.2. edgeR analysis was carried out in R version 3.5.1. A volcano plot was created to visualize DEGs. Custom scripts are available on GitHub at link https://github.com/mclemente-igenomix/garrido_et_al_2021.

Transcriptomic fingerprinting definition and validation

Genes with assigned EntrezID with an FDR cutoff of 0.05 and an expression ≥1.4-fold higher in the sPE vs. control training set samples were selected to define a fingerprint associated with DD in sPE. Targeted analysis of fingerprinting genes was performed using the validation set of samples. PCA and unsupervised hierarchical clustering with a Canberra distance based on gene signature were performed comparing sPE to control specimens. Custom scripts are available on GitHub at https://github.com/mclemente-igenomix/garrido_et_al_2021.

Enrichment analysis

Gene Ontology (GO) analyses were conducted to obtain biological processes using the goana function in edgeR (Robinson et al., 2010). The input genes were those 120 included in the fingerprinting (Figure 3—source data 1). The p-value adjustment method was FDR with a cutoff of 0.05 (FDR < 0.05; Figure 3—source data 2). Custom scripts are available on GitHub at https://github.com/mclemente-igenomix/garrido_et_al_2021.

Interaction network

An interaction network between proteins encoded by DD fingerprinting genes was created using the functional analysis suite String (Jensen et al., 2009). To construct the network, the interactions included were from curated databases and included experimentally determined and predicted interactions, textmining, co-expression information. The clustering algorithm k-means was applied based on the distance matrix obtained from the String global scores. The network was visualized using Cytoscape software (Shannon et al., 2003). Hub genes were extracted using the maximal clique centrality (MCC) and maximum neighborhood component (MNC) of the cytoHubba plugin (Chin et al., 2014). The overlapping genes identified by the two topological analysis methods were selected as the hub genes.

qRT-PCR

Gene expression of IHH, MSX2, ESR1, and PGR isoforms in the endometrial tissue from a subset of women with prior sPE (n = 13) compared to controls (n = 9) was obtained by RT-qPCR. Specific primers for each gene are described in Supplementary file 3. cDNA was generated from 400 ng of RNA using the SuperScript VILO cDNA Synthesis Kit (Thermo Fisher Scientific). Template cDNA was diluted 5 in 20 and 1 µL was used in each PCR. Real-time PCR was performed in duplicate in 10 µL using commercially validated Kapa SYBR fast qPCR kit (Kapa Biosystems Inc, Basilea, Switzerland) and the Lightcycler 480 (Roche Molecular Systems, Inc, Pleasanton, CA) detection system. Samples were run in duplicate along with appropriate controls (i.e., no template, no RT). Cycling conditions were as follows: 95°C for 3 min, 40 cycles of 95°C for 10 s, 60°C for 20 s, and 72°C for 1 s. A melting curve was done following the product specifications. Data were analyzed using the comparative Ct method (2−∆∆CT). Data were normalized to the housekeeping gene β-actin, and changes in gene expression were calculated using the ΔΔCT method with the control group used as the calibrator; values are illustrated relative to median in the control group. The relative expression of PGR-A mRNA was calculated by subtracting the relative expression of PGR-B mRNA from that of PGR total.

Immunofluorescence of tissue sections

Endometrial tissue samples were fixed in 4% paraformaldehyde and preserved in paraffin-embedded blocks. For immunostaining, tissue sections were deparaffinated and rehydrated. Antigen retrieval was performed with buffer citrate 1× at 100°C for 10 min. Then, non-specific reactivity was blocked by incubation in 5% BSA/0.1% PBS-Tween 20 at room temperature for 30 min. Sections were incubated at room temperature for 1.5 hr with primary antibodies (1:50 rabbit monoclonal anti-human progesterone receptor, Abcam, Cambridge, UK) and 1:50 mouse monoclonal anti-human estrogen receptor 1 (Santa Cruz Biotechnology, CA) diluted in 3% BSA/0.1% PBS-Tween 20. Then, slides were washed two times for 10 min with 0.1% PBS-Tween 20 before they were incubated for 1 hr at room temperature with AlexaFluor-conjugated secondary antibodies diluted in 3% BSA/0.1% PBS-Tween 20 (1:1000). Finally, slides were washed two times in 0.1% PBS-Tween 20. To visualize nuclei, 4′,6-diamidino-2-phenylindole at 400 ng/µL was used. Tissue sections were examined using a EVOS M5000 microscope.

Statistical analysis

Clinical data are expressed as mean ± standard error mean (SEM). Clinical data were evaluated by Wilcoxon test for comparisons between sPE and control samples. Statistical significance was set at p<0.05. Differential expression analysis was performed using the R package edgeR.

Results

Endometrial transcriptome alterations during decidualization in sPE

To identify transcriptomic alterations during decidualization in sPE, we applied global RNA sequencing (RNA-seq) to endometrial biopsies obtained in the late secretory phase from women who developed sPE in a previous pregnancy (n = 24) and controls who never had sPE (n = 16) (GSE172381). Clinical maternal and neonatal characteristics of the participants are summarized in Supplementary file 1. After quality trimming and filtering, reads were aligned to the reference genome hg19. The 40 samples produced 56,638 raw sequencing genes; after normalization, 18,301 genes were included in the analysis. Biological and technical variables for each donor were considered to discard confounding effects on the transcriptomic profile (Supplementary file 2). Controls included women who had a preterm birth with no signs of infection (n = 8) and women who gave birth at full term with normal obstetric outcomes (n = 8). Transcriptomic profiles were compared by differential expression analysis, revealing no significant changes in the endometrial transcriptome between preterm and term controls (FDR ≥ 0.05; Figure 1—figure supplement 1A). Principal component analysis (PCA) supported that there was no underlying pattern of distribution depending on gestational age at delivery (Figure 1—figure supplement 1B). Once we ruled out bias on controls, we randomly split samples into two cohorts, a training set (70%) and a test set (30%) (Figure 1A). Random sampling occurred within each class (sPE and controls), so overall class distribution of the data was preserved. The training set (n = 29) was used for the identification of molecular fingerprinting encoding DD in sPE, while the test set (n = 11) was used to confirm our findings. All samples in both cohorts were processed and sequenced in the same manner.

Figure 1 with 1 supplement see all
Global RNA-seq transcriptomic results revealed 593 differentially expressed genes (DEGs) in severe preeclampsia (sPE) vs. control samples.

(A) Schematic drawing of the study design used to identify and validate defective decidualization (DD) fingerprinting in sPE. (B) Statistical significance (-log10 FDR) vs. gene expression log2 fold change (FC) is displayed as a volcano plot of global RNA-seq results. Label indicates: downregulated in sPE (blue dots); upregulated in sPE (red dots); not significant genes (grey dots). (C) Heatmap showing the 25 most upregulated and downregulated genes (total = 593; Figure 1—source data 1) of control vs. sPE samples. See also Figure 1—source data 1.

Figure 1—source data 1

The 593 statistically differentially expressed genes (false discovery rate [FDR] < 0.05) with at least 1.2-fold change (FC ≥ 1.2) in severe preeclampsia (sPE) vs. control cases obtained from RNA-seq analysis.

https://cdn.elifesciences.org/articles/70753/elife-70753-fig1-data1-v1.xlsx

Transcriptional analysis in the training set was performed by comparing gene expression patterns in sPE (n = 17) and controls (n = 12). This comparison revealed 593 DEGs based on FDR < 0.05 and with at least 1.2 FC between groups (FC ≥ 1.2). DEGs are shown in the volcano plot through yellow dots (Figure 1B). A total of 155 upregulated and 438 downregulated DEGs were identified as being associated with DD in sPE (Figure 1C; complete list in Figure 1—source data 1). Downregulated transcripts include those involved in decidualization, such as MMP3, PRL, IL-6, and IHH; and genes associated with signaling (e.g., NR4A3 and IL8), growth factors (e.g., FGF1 and FGF7), angiogenesis (e.g., EDN2 and TMEM215), and immune response (CCL20, CXCL3, and IGHG1). Upregulated genes are involved in amino acid metabolic/catabolic processes (IDO2 and CAPN3), transport, and oxidoreductase activity.

Comparison of DD transcriptomics in previous sPE in vivo vs. in vitro

We previously described DD in human endometrial stromal cells (hESCs) isolated from women with previous sPE compared to women with normal obstetric outcomes, but this finding was restricted to the stromal cell population using an in vitro decidualization cell culture model (Garrido-Gomez et al., 2017). Here, we compared DD overlapping between DEGs reported in vitro (n = 129) vs. in vivo (n = 593) in sPE compared to control women. Nine genes were overlapped between the two datasets (Figure 2A); one gene was upregulated (ERP27), and eight genes were downregulated (e.g., ISM1, MEST, MFAP2, and REEP2). The expression pattern of common genes is presented as a box plot using counts per million, corroborating significant differential expression between sPE and control (Figure 2B). Recently, in vivo transcriptomics of endometrium at single-cell resolution across the menstrual cycle were characterized (Wang et al., 2020). Transcriptome profiles of stromal fibroblasts from the late secretory phase allowed the identification of deregulated genes in sPE as associated to hESC. We found that 263 genes from the 593 DEGs in sPE vs. control are expressed by hESC (Figure 2C). Taken together, the in vivo assessment provides a broad spectrum of dysregulated transcripts comparing with previous in vitro findings, which includes a high concordance with in vivo hESC genes.

Defective decidualization (DD) transcriptomics in vitro vs. in vivo.

(A) Common genes between previous in vitro (left) and current in vivo approaches analyzing decidualization (right). Nine genes overlap in both approaches. (B) Box plot showing the average expression of the nine common genes between control (blue boxes) and severe preeclampsia (sPE) (orange boxes) samples. (C) From the 593 differentially expressed genes (DEGs) obtained by global RNA-seq, a subset of 263 DEGs were identified as genes with a human endometrial stromal cell (hESC) origin using the scRNA-seq data published by Wang et al., 2020. See also Figure 2—source data 1.

Figure 2—source data 1

The 593 statistically differentially expressed genes (false discovery rate [FDR] < 0.05) with at least 1.2-fold-change (FC ≥ 1.2) in severe preeclampsia (sPE) vs. control cases obtained from RNA-seq analysis.

https://cdn.elifesciences.org/articles/70753/elife-70753-fig2-data1-v1.xlsx

Identification of the fingerprint encoding human endometrial DD

To formulate the transcriptomic signature that encodes DD detected in sPE in vivo, we selected genes with significant dysregulation (FDR < 0.05) and at least 1.4-fold increase (FC ≥1.4) between sPE and control with assigned EntrezID. A volcano plot shows 120 DEGs meeting these criteria included in the final DD signature (Figure 3A; complete list of genes is included in Figure 3—source data 1).

Severe preeclampsia defective decidualization (sPE-DD) fingerprint composed of 120 differentially expressed genes (DEGs).

(A) Volcano plot showing downregulated (blue) and upregulated (red) genes in sPE from the DD fingerprint. Each point represents one gene; gray points are the rest of the genes obtained in the global RNA-seq analysis. (B) The three most highly downregulated biological process for each major category (red, cell cycle; yellow, DNA damage response; green, cell signaling; blue, cellular response; gray, cell motility; purple, extracellular matrix; pink, immune response; brown, reproductive process). Enrichment index was calculated by -log(p-value). (C) Clustering of DD fingerprint genes shown for reproductive process, response to bacterial molecules, extracellular matrix organization, regulation of receptor signaling, and response to hormones. See also Figure 3—source data 1 and Figure 3—source data 2.

Figure 3—source data 1

List of genes selected as defective decidualization signature in severe preeclampsia (sPE) (120 differentially expressed genes [DEGs] with false discovery rate [FDR] < 0.05 and fold-change [FC] ≥ 1.4).

https://cdn.elifesciences.org/articles/70753/elife-70753-fig3-data1-v1.xlsx
Figure 3—source data 2

Biological process Gene Ontology (GO) terms computed by the 120 genes included in the defective decidualization (DD) fingerprinting in severe preeclampsia (sPE) (N, number of genes associated to GO term; DE, number of genes differentially expressed in this GO term; P.DE, unadjusted p-value; FDR, adjusted p-value).

https://cdn.elifesciences.org/articles/70753/elife-70753-fig3-data2-v1.xlsx

GO analysis of the gene signature associated with DD in sPE identified 151 enriched biological processes downregulated (FDR < 0.05). These pathways were associated with cell cycle, DNA damage response, cell signaling, cellular response, cell motility, extracellular matrix, immune response, and reproductive process (Figure 3B). All are hallmarks of impaired decidualization and sPE pathogenesis. We identified fingerprinting genes representative of the altered pathways in sPE, such as IL6 and TNF, regulating the response to bacterial molecules, MMP3 and MMP1 participating in the extracellular matrix organization, and TNF, IL8, and FGF1 implicated in the downregulated receptor signaling (Figure 3C). Functional analysis revealed that the 120 DEGs included in DD fingerprinting are implicated in pathways related to decidualization, corroborating the maternal contribution to sPE. Interestingly, the number of downregulated genes was higher than the number of upregulated genes in sPE compared to controls, suggesting that, in vivo, DD may be induced by the lack of expression of a subset of genes.

Based on the 120 genes included in the DD signature, PCA showed that sPE and control samples clustered separately in two groups, except for three control samples (C20, C21, and C22) (Figure 4A). High variance between groups was effectively captured in the first two principal components. Unsupervised hierarchical clustering analysis confirmed that gene fingerprinting effectively segregated the two groups: one encompassing mainly controls and the other mainly sPE samples (Figure 4B). The same three controls clustered with the sPE group, recreating the PCA results.

Validation of the defective decidualization (DD) fingerprint in severe preeclampsia (sPE).

(A) Principal component analysis (PCA) based on 120 genes included in the fingerprinting in the training set. Each sample is represented as a colored point (blue, control; orange, sPE). (B) Heatmap dendrogram of expression of the 120 genes included in the final fingerprinting for each sample of the training set (control, n = 12; sPE, n = 17). (C) PCA based on the fingerprinting in the test set. Each sample is represented as a colored point (blue, control; orange, sPE). (D) Heatmap dendrogram of expression of the 120 genes included in the final fingerprinting for each sample of the test set (control, n = 4; sPE, n = 7). See also Figure 4—source data 1.

Figure 4—source data 1

List of genes selected as defective decidualization signature in severe preeclampsia (sPE) (120 differentially expressed genes [DEGs] with false discovery rate [FDR] < 0.05 and fold-change [FC] ≥ 1.4).

https://cdn.elifesciences.org/articles/70753/elife-70753-fig4-data1-v1.xlsx

To validate the DD gene signature in an independent cohort of samples (sPE [n = 7] vs. control [n = 4]), PCA based on these transcripts effectively segregated samples in two homogeneous groups (Figure 4C), corroborated by hierarchical clustering (Figure 4D). These genes successfully grouped 100% of controls and 85.7% of sPE cases supporting DD in sPE.

DD fingerprint in sPE is connected to ER1 and PR-B

Of the 120 genes in the DD signature, 94 endometrial enriched genes encode for specific proteins reported by the Human Protein Atlas (Uhlén et al., 2015). Interestingly, 45 of those genes (47.9%) were included in the transcriptome modulated by ESR1 (Hewitt et al., 2010), and 43 genes (45.7%) overlapped with the transcriptome and cistrome associated with PGR (Mazur et al., 2015; Figure 5A). Regarding target genes of ER1 and PR, the database of Human Transcription Factor Targets (hTFtarget) reported 17 genes responsive to ER1 and 50 target genes modulated by PR, based on epigenomic, CHIP-seq, or motif evidence (Zhang et al., 2020).

Estrogen receptor 1 (ER1) and progesterone receptor-B (PR-B) are linked to defective decidualization (DD) fingerprinting in severe preeclampsia (sPE).

(A) Venn diagram displaying genes included in the fingerprinting (120) predominantly expressed in the endometrium based on Human Protein Atlas data that overlap with genes modulated by ESR1 described by Hewitt et al., 2010 and genes associated with PGR silencing described by Mazur et al., 2015. (B) Network showing the connections between proteins codified by DD fingerprinting and the hormonal receptors, ER1 and PR. Shapes indicate different clusters established by String k-means method. Squares, cluster involved in gland morphogenesis and cell migration; circles, cluster involved in extracellular matrix organization and stromal cell differentiation; hexagons; cluster involved in cellular response to DNA damage and regulation of cell cycle. Color gradient indicate gene expression in terms of log2FC. Hub genes are shown with an asterisk. (C-H) Gene expression levels of IHH, MSX2, ESR1, PGR, PGR-A, and PGR-B assessed for sPE (n=13) vs. controls (n=9) by RT-qPCR (gray bars, control; green bars, sPE). RT-qPCR values are expressed as mean± SE. *** p<0.001, ** p<0.01, *p<0.05. (I-J) Tissue sections of control (n=4) and sPE (n=4) endometrium during late secretory phase were immunostained with antibody against ER1 or PR. Nuclei were visualized with DAPI. Scale bar: 50 µM.

We evaluated the interaction between steroid receptor signaling and the proteins encoded by DD fingerprinting genes in sPE by building a dynamic network including ER1 and PR. String software (Jensen et al., 2009) was used to construct network connections visualized with Cytoscape software (Shannon et al., 2003). The interactome contained 117 nodes directly interconnected by 361 edges (Figure 5B). This DD fingerprint network showed a highly enriched protein–protein interaction (PPI) in sPE; indeed, the interconnection between nodes was significantly higher than the 93 edges expected (PPI enrichment p<1.0e-16). Clustering revealed three main modules based on their connectivity degree, with functionally relevant genes involved in gland morphogenesis, cell migration, extracellular matrix organization, stromal cell differentiation, cellular response to DNA damage stimulus, and regulation of cell cycle. The hub genes were determined by overlapping the top 10 genes obtained using two topological analysis methods in the cytoHubba plugin (Chin et al., 2014), MCC, and MNC. Five genes were selected, all of which were downregulated. Interestingly, both ER1 and PR were strongly embedded in the network and highly connected with DD fingerprinting, highlighting the interaction of hormonal receptors with notable decidualization mediators such as IHH and MSX2 validated by RT-qPCR (Figure 5C and D). Furthermore, the interactome demonstrated a direct interaction between ER1 and PR. These results support the transcriptomic dysfunction of the genes present in the DD signature through imbalanced hormone receptor signaling in sPE.

We then analyzed the expression of ESR1 and PGR in the endometrial tissue from a subset of women with prior sPE (N = 13) compared to controls (N = 9) by RT-qPCR. We found reduced expression of transcripts encoding the hormone receptors ESR1 (p<0.05) and PGR (p<0.001) in sPE patients (Figure 5E and F). In-depth expression analyses revealed that the isoform PGR-B was significantly downregulated in sPE vs. controls (p<0.01), while the isoform PGR-A was unaffected (p>0.05) (Figure 5G and H). These results were confirmed at the protein level by immunohistochemical analysis of ER1 and PR-B in endometrial biopsies collected in the late secretory phase from women with previous sPE (n = 4) and controls (n = 4) (Figure 5I and J). Both receptors were highly expressed through the decidualized endometrium, especially in the secretory glands in controls. In contrast, their expression was greatly reduced or absent in sPE samples. These results suggest that the DD transcriptomic signature implicates dysregulated ER1 and PR-B signaling in the late secretory phase in sPE patients.

Discussion

Most scientific and clinical diagnostic efforts in sPE focus on placental surrogates. In this context, shallow cytotrophoblast invasion induces deficient vascular remodeling and ultimately aberrant placentation. This leads to placental ischemia and the release of soluble factors that induce the maternal syndrome, including the imbalanced levels of soluble fms-like tyrosine kinase 1 (sFLT1) and placental growth factor (PlGF) (Rana et al., 2019; Powe et al., 2011). sFLT1 protein binds to PlGF, preventing its interaction with endothelial receptors and leading to endothelial dysfunction. Accordingly, sFLT1 is increased, while free PlGF is decreased in serum from women with PE (Levine et al., 2004). The sFLT1/PlGF ratio has been proposed as a biomarker showing a positive predictive value of 36.7%, with 66.2% sensitivity and 83.1% specificity but its application is effective only 4 weeks before PE symptoms manifest (Zeisler et al., 2016). Therefore, during the first trimester there is a lack of high sensitivity and specificity screening methods to detect sPE early and prevent mortality and morbidity.

Current strategies based on placental dysfunction provide delayed results for preventive interventions and new approaches are urgently needed. Recent studies suggest that sPE might also be a disorder of the decidua, opening new avenues. In this sense, there may be a molecular signature associated with impaired endometrial maturation during early pregnancy in women who develop sPE (Rabaglino et al., 2015). Our previous work demonstrated that hESC isolated from patients with previous sPE failed to decidualize in vitro, suggesting a role of the maternal factor in the development of this disease (Garrido-Gomez et al., 2017). In addition, molecular pathways of dysregulated decidualization in PE are also found in endometrial disorders such as implantation failure, recurrent miscarriage, and endometriosis (Rabaglino and Conrad, 2019). Increasing evidence supports the detriment of inappropriate decidualization before pregnancy to reproductive outcomes, and the role of DD in the origin of sPE (Garrido-Gómez et al., 2020; Ng et al., 2020).

In the present study, we highlighted the underlying molecular defect that may explain in vivo decidualization failure as an important contributor to shallow placental invasion in sPE. For this purpose, we investigated the role of the decidua by leveraging global RNA-seq in late secretory endometrium. We identified 593 genes differentially expressed in sPE compared to controls, including genes involved in decidualization such as PRL, IL6, and IHH, as well as several novel genes. Nine of the DEGs overlapped with DEGs obtained in our previous in vitro decidualization study. In vivo endometrial biopsies include other cells in addition to hESC, and thus the expected degree of overlap between both in vitro and in vivo approaches should be modest as such was observed. However, we identified a large percentage of those DEGs that were associated with in vivo transcriptomic profile of hESC resolved at cell level from a healthy late secretory endometrium. Thus, the high number of our identified in vivo DEGs reflects the high complexity of decidualization at cellular heterogeneity in a physiological stage of women with previous sPE.

Previous reports of the decidual transcriptome in PE and sPE (Garrido-Gomez et al., 2017; Løset et al., 2011) revealed the gene expression profile associated with the condition at the time of delivery. Here, we analyzed samples collected years after the affected pregnancy; thus, it could be interesting to find dysregulated genes in common among these approaches. We compared our transcriptomic results and the previously reported dysregulated decidual genes, obtaining 1.3% (Garrido-Gomez et al., 2017) and 1.7% (Løset et al., 2011) overlap in affected genes. These discordances may not be unexpected and are consistent with recent results (Rabaglino and Conrad, 2019). The decidua basalis transcriptome at delivery was compared with the transcriptome of decidua at ⁓11.5 gestational weeks and in vitro decidualized hESC from women who experienced sPE. Both analyses revealed little, if any, overlap between molecular signatures. In contrast, the signature encoding in vitro DD of hESC years after pregnancy overlapped significantly with the dysregulated profile found in decidual samples at the beginning of pregnancy (Rabaglino and Conrad, 2019). Thus, decidual gene expression patterns during the clinically active disease largely differ from those observed during endometrial decidualization at the end of the menstrual cycle and early pregnancy, perhaps reflecting a consequence rather than the origin of sPE.

From the 593 genes differentially expressed in sPE compared to controls, we identified a DD signature comprising 120 genes associated with sPE. This sPE-DD signature includes genes that allowed us to segregate samples from the training set into sPE and control groups, which were confirmed by the test set. One sPE was misclustered in the test set, such that 90.9% of samples from an independent cohort were properly clustered in the dendrogram. Having controlled for confounding effects of biological and technical variables, we consider that this misclustering is consistent with the nature of decidualization and its inherent variability. Decidualization is a highly dynamic process governed by (1) inter-individual variability in the endometrial menstrual cycle supported by displacement of the window of implantation in one out of four patients experiencing recurrent implantation failure (Wang et al., 2020; Diaz-Gimeno et al., 2011), (2) the physiology of the spatial expansion of decidualization process starting in some areas around spiral arteries and extending to the entire endometrium during the last days of the menstrual cycle (Gellersen and Brosens, 2014), (3) and the random spatial sampling during the endometrial biopsy that could influences cell-type proportions, which is inherent to the experimental strategy used in our investigation. Decidua sample segregation is consistent with the results presented by Munchel et al., 2020, who classified PE patients based on circulating RNA (C-RNA). In our study, we determined that the gene expression associated with DD had the highest potential to segregate sPE.

Interestingly, most of the DD fingerprint genes in sPE were related to the downregulation of ESR1 and PGR, specifically PGR-B. We hypothesize that low expression of PGR-B and ESR1 activates endometrial decidualization by dysregulating progesterone (P4) and estrogen (E2) action. Consequently, P4- and E2-related cellular signaling may be compromised, leading to the development of DD. Stromal cell differentiation, stromal–epithelial crosstalk (Wang et al., 2017), extracellular matrix degradation (Itoh et al., 2012), immune system response, and endothelial function (Okada et al., 2018) may all be disrupted by low expression of PGR-B and ESR1. Remarkably, the local immune system appeared to be dysregulated due to the altered expression of interleukins, cytokines, chemokines, and immunoglobulin, which could disrupt the tolerant microenvironment at the maternal–fetal interface in pregnancy and contribute to the development of sPE (Erlebacher, 2013; Harris et al., 2019).

Based on our findings, we postulate that in non-sPE pregnancies balanced hormonal signaling leads to proper decidualization, which in turn interacts with immune and endothelial cells to control cytotrophoblast invasion (Figure 6A). P4 and E2 activate their receptors in the epithelium and signal to the stromal compartment. Likewise, stromal PR is activated by ER1, which induces target genes involved in decidualization such as MSX2, CCNB1, and IL6. In contrast, in sPE, endometrial ESR1 and IHH are downregulated, decreasing PGR expression (Figure 6B) and leading to compromised decidualization, endothelial dysfunction, and local immune dysregulation.

Modeling of the molecular mechanism for defective decidualization (DD) in severe preeclampsia (sPE).

(A) Decidualization induced by P4 and E2 in control pregnancy including the interaction of immune response and endothelium. (B) Hypothetical network that could link DD and dysregulated hormone signaling in sPE. All genes were downregulated. Biological processes specified are candidates to be impaired based on functions associated with the observed dysregulation. Red arrows show the downregulation of decidualization modulators.

In sPE, we found evidence of impaired epithelial–stromal crosstalk through IHH and FGF1, which could lead to an imbalance of hESC proliferation and differentiation (Cha et al., 2012; Wang et al., 2018). In support of this, we observed downregulated genes involved in stromal cell differentiation including IL6, CCNB1, CCNB2, and MSX2. Further, metalloproteinases such as MMP1, MMP7, and MMP11 were downregulated in sPE, which could be associated to defective cell motility and extracellular matrix degradation. Other transcriptomic deregulations in sPE are connected to endothelial and immune system dysfunction, such as EGR3 and MMP3, both of which are associated with altered angiogenesis (Chen et al., 2020; Frieling et al., 2020); TNF is involved in endothelial inflammation disbalance (Marcos-Ramiro et al., 2014) and influences the tolerant microenvironment at the maternal–fetal interface (Li et al., 2018). Taken together, our findings reveal a DD transcriptomic fingerprint in sPE driven by an imbalance in ER1 and PR signaling.

Our findings reveal significant gene expression dysregulation underlying DD in the late secretory phase in women who have had sPE. The potential origin of sPE may lie in the downregulation of ESR1 and PGR-B. Both receptors are strongly coordinated to regulate decidualization and PGR expression is induced by ER1, which is inhibited by PR (Patel et al., 2015). E2 and P4 act through ER1 and PR in the epithelium and stroma and modulate the transcriptome and to promote crosstalk between both compartments (Winuthayanon et al., 2010). ER1 in the epithelium regulates stromal decidualization via paracrine mechanisms mediated by leukemia inhibitory factor (LIF), which controls IHH expression that transduce the signal activation of PR in the stroma (Pawar et al., 2015). Also, stromal ER1 activates PR in the same compartment (Kaya Okur et al., 2016). These findings could inform the development of therapeutic targets to restore optimal decidualization in sPE. Further studies are needed to fully elucidate the hormone signaling pathways that become dysregulated in sPE and the role of DD in the manifestation of the condition.

Preeclampsia is a syndrome in which different condition types might coexist; here, we focused on the maternal contribution to sPE through DD. We are not claiming that our findings are related to sPE heritability (DNA), but rather to the pathogenesis of the decidualization defect (RNA expression) involved in the origin of PE (Garrido-Gómez et al., 2020; Ng et al., 2020). Our findings reinforce a maternal cause for sPE through DD. This condition may result from an aberrant response to estrogen and progesterone mediated by ER1 and PR-B signaling. However, the primary driver of the predisposition to undergo decidualization resistance and its link with the main risk factors of sPE remain to be determined. Our work is an important step toward the development of new strategies that enable early assessment of risk for sPE and might prompt new therapeutic strategies to treat this enigmatic pathological condition. Future studies should focus on the translational potential of the DD fingerprinting to develop new noninvasive strategies based on circulating RNA to improve diagnosis and prognostication for women with sPE.

Data availability

Transcriptomic data were deposited in the Gene Expression Omnibus database (accession number GSE172381) (Included in results and material and methods section). Custom scripts are available on GitHub at link https://github.com/mclemente-igenomix/garrido_et_al_2021, copy archived at https://archive.softwareheritage.org/swh:1:rev:c5a81d119229bc9d46c127c039a85de327f9b9c5.

The following data sets were generated
The following previously published data sets were used
    1. Wang W
    2. Vilella F
    3. Alama P
    4. Moreno I
    5. Mignardi M
    6. Isakova A
    7. Wenying P
    8. Simon C
    9. Quake SR
    (2020) NCBI Gene Expression Omnibus
    ID GSE111976. Single cell RNA-seq analysis on human endometrium across the natural menstrual cycle.

References

Decision letter

  1. Bérénice A Benayoun
    Reviewing Editor; University of Southern California, United States
  2. Kathryn Song Eng Cheah
    Senior Editor; The University of Hong Kong, Hong Kong

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Acceptance summary:

Garrido-Gomez et al., generate a unique dataset profiling the transcriptomes late-stage endometrium of women with prior pregnancies leading to severe preeclampsia (sPE) compared to that of non-preeclamptic pregnancies. Although the question of molecular drivers preeclampsia has been explored at the molecular level directly in the placenta, this study provides a unique view in the endometrium of women years after the pregnancy in question, which is consistent with the fact that a previous preeclamptic pregnancy is one of the best predictors for a subsequent preeclamptic pregnancy.

Decision letter after peer review:

Thank you for submitting your article "Disrupted PGR-B and ESR1 signaling underlies preconceptional defective decidualization linked to severe preeclampsia" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Kathryn Cheah as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

Although the 3 reviewers noted the potential of the study, major concerns were raised about the soundness/clarity of the transcriptomics analysis, as well as lack of clarity on the study population. Specifically:

1. The manuscript needs to undergo major revisions in the statistics and transcriptomics analyses (see specific points with details in reviewers 1 and 2 comments):

a. The use of fold-change thresholding needs to be removed/updated.

b. Multiple hypothesis correction needs to be used prior to determining significance with genome-wide datasets. Please use FDR correction or other accepted methods.

c. Only use Student t-tests when data normality can be confirmed, otherwise, use non-parametric tests.

d. Since this study is done on endometria sampled years after actual pregnancy, it would be crucial to compare to public datasets on actual placentas, to determine overlap and differences.

e. PCA cannot be used to determine clustering or lack thereof, as well as whether differences exist or not.

f. There are concerns about batches and co-variates not taken into account (see reviewers 1 and 3 points) which could lead to the absence of DEGs. These need to be clearly stated and provided in a table. They also need to be corrected for, either using batch removal and/or multivariate modeling.

g. A continuous GSEA-type analysis would be more informative than a thresholded approach for functional enrichment (i.e. GO) analysis.

h. RT-qPCR (especially on the exact same sample set) is not an appropriate validation for robust RNA-seq datasets.

2. The authors needs to substantially revise and improve the details on methodology (i.e. RIN threshold), as well as deposit all code and data to public databases for re-analysis. This is crucial for long-term reproducibility and use of this as a resource. Analytical choices against the general rule of the field (such as using hg19 and not the more recent better annotated hg38) should be justified or at least explicitly addressed.

3. The patient population needs to be better defined (see reviewer 3 points, and to a lesser degree reviewer 1).

4. Make sure to define all abbreviations, and avoid using non-standard abbreviations to improve readability.

Reviewer #1 (Recommendations for the authors):

1. There are major concerns about the statistics in the analysis of the dataset.

a. Fold-change thresholding is used in multiple places in the analysis (e.g. line 94, 112, 128, etc.). This is a big problem, as it is known to lead to poor FDR control in the case of statistical tests not designed to take fold change into account (like edgeR used here). Thus, the practice of fold-change filtering without FDR control is considered problematic by statisticians (see PMID: 19176553). Thus, either the authors need to update their results after removing this problematic criteria to properly control for FDR, or all analyses should be rerun with a method that takes fold-change into account for FDR control (e.g. the TREAT method PMID: 19176553).

b. For the GO enrichment analysis (line 133), the authors only indicate a p-value filter (p < 0.005). Was there not implementation of multiple testing correction? The data should be re-analyzed after proper multiple hypothesis correction control (native DEseq2 FDR, Bonferroni, or other). In addition, what was the gene list background used to compute enrichment? The nature of the background has been shown to hugely impact results and needs to be explicitly detailed in the methods.

c. The authors reports using Student t-tests (section starting line 404). Unless normality is proven (e.g. with a Shapiro-Wilkes test), parametric tests that relie on normal distributions should never be used. We recommend that all tests be rerun with non-parametric equivalents (e.g. Wilcoxon or Mann Whitney) to guarantee these assumptions are not violated.

2. With human sample studies, it is unlikely that all samples were processed exactly in parallel (i.e. all samples collected the same day, kept in stabilizer for the same amount of time, etc.), which all but guarantees the existence of batch effects. There are also a number of stated parameters (i.e. age of the donor, time since pregnancy 1-8 years, day in cycle 22-32, etc.) that could influence results without a relevant biological driver.

a. In addition to an aggregated table (Table 1), the authors should generate a table for each anonymized sample for key characteristics so that the resource can be analyzed with these caveats considered both in this current paper and to facilitate all future reanalysis attempts. This should include for each sample individually:

i. age of the donor.

ii. time since last pregnancy.

iii. day of cycle for biopsy.

iv. date of processing/batch (including biopsy batch, library batch AND sequencing batch).

v. RIN of RNA sample.

vi. Sequencing depth (raw reads per library).

vii. Percent reads mapping to the genome reference (here hg19).

viii. Ethnicity of donor.

b. The authors should implement a multivariate model to account for all these potential confounding effects so as to focus on the core biological signal. In addition, the use of a software like SVA or RUVseq is recommended to remove unwanted technical variation prior to analyses.

3. Additional methodological information is needed for long-term reproducibility of analyses.

a. For reproducibility of code and analyses, all analytical code for this study should be deposited in a repository such as Github or made available as a Supplementary file.

b. Please include all version numbers for all used software (e.g. R, etc.) packages and R packages (e.g. edgeR, goana, etc.), as well as all command parameters where relevant. The same should be applied to annotation databases (i.e. GO used here), and if a version number doesn't exist, date of access should be provided.

c. There is no mention of the softwares (nor versions) used to trim RNA-seq reads, to map to the hg19 genome or to aggregate counts to genes. All used software should be clearly named and their use described with all necessary options for reproducibility.

d. Line 84-88/Figure 1A, were the subtypes of "control" samples (pre vs. full-term) equally distributed in training vs. testing? Indeed, although the authors argue that the subtype made no difference on transcriptional profiles (i.e. Sup Figure 1A PCA), this reviews sees a general (if not clean separation) on PC1, that would probably be enhanced if sources of spurious variation (see point #2) were taken into account. In any case, the data provided cannot be used to support the statement "demonstrated an absence of clustering" (line 80), since the statistical tests used are not meant to prove the null hypothesis but to reject it when possible (i.e. confusion between type I and II error). Please amend the analysis to reflect this.

e. cDNA library cannot be validated for RIN (line 330), as only total RNA samples can be. The entire method section needs to be completely rechecked and rewritten for accuracy to avoid this kind of errors.

f. All used QC filters should be stated explicitly in the methods (see lines 338-342).

4. Since RT-qpCR is known to be less sensitive and more prone to normalizing biases (due to the choice of control genes which may vary themselves) than RNA-seq (which performs unbiased normalization at the transcriptome-wide level), it is unclear why validation was performed on the same samples that were processed by RNA-seq (line 375-389). Unless the authors want to include RT-qPCR data on an INDEPENDENT cohort of patients, these results are circular and shouldn't be included in a revised manuscript. If RT-qPCR on a new cohort is included, make sure to include information on which normalizing amplicons were used for δ CT calculation.

5. Can the authors explain why 3 specific samples (C20, 21, 22; line 147-152) do NOT cluster appropriately with their "signature"? Did they have any biological specificities (see need for sample-by-sample information as raised in point #2)? The fact that 3 control samples cluster with the sPE samples somewhat invalidate the signature as a signature of sPE, and may be the result of improper correction for batching or other technical (or irrelevant biological) noise. Please discuss explicitly what may be happening here in a revised manuscript.

Reviewer #2 (Recommendations for the authors):

In sum, I feel that this paper brings some interesting insights on the decidua status and preeclampsia; the degree of novelty is nevertheless unclear to me, and the clinical application seems far-fetched. There are several points that may be better presented and discussed.

Detail of the recommendations to the authors:

1. The study of the decidual transcriptome in preeclampsia has also already been performed by other teams in the past, rather at the moment of the disease than later. For instance, the paper of Mari Loset (Am J Obst Gynecol, 2011), is, I feel, an important base for comparison with the results presented here. So much so, that it used the same type of expressional approaches, and detected differential transcripts and cascades. I feel that the authors should analyze their own datasets in light of this type of seminal papers and discuss the commonalities and differences found.

2. In the overrepresented pathways in the Loset's work, the regulation by ESR1 and PGR was not obvious. Does it mean that analyzing samples 3-4 years after the event (normal or pathological pregnancy) leads to remnants of the disease at the uterine level, or is there a genetic predisposition? This is a question that cannot be truly addressed by the present dataset, and is important if the objective is to define markers, as stated below.

The introduction gives the necessary details to understand the question raised, finding a signature that could help assessing the risk of preeclampsia.

3. However, the genetic part of preeclampsia is estimated by a heritability of ~55%. The genetic decomposition of this heritability presented by Cnattingius and coworkers (2004) indicate that only part of this genetics is connected to the maternal genetic background (which concerns two organs: the uterus and thus the decidua) and the placenta for half of its genetics. Part of the genetics of preeclampsia is associated only to the paternal genetics (estimated at ~1/3 of the heritability), therefore, it is clear that the risk cannot systematically be found by analysis of the female (uterine-decidua) expression profile only. Alternatively, the difference observed by the authors may be a consequence of the previous preeclampsia not based on a predisposition, but if this is the case, it limits considerably the usefulness of the markers found since preeclampsia, and in particular severe preeclampsia is at 75% a disease of the first pregnancy. The authors should recognize this as an important limit as an early prenatal screening strategy.

The initial analysis is based upon the endometrial tissue of 24 women that had a severe preeclampsia, 16 controls (8 preterm and 8 term births). Gestational age at delivery had apparently no influence on the 'control' groups. Since the samples were not collected at the end of gestation, which means that the occurrence of preterm may be connected to the placenta and not to the maternal uterine situation.

4. Nevertheless, there are only 8+8 samples so it is not possible to prove that in some case a decidual expression is not involved. In addition, using PCA to check that there is no clustering is not enough to certify that there are no differentially expressed genes that separate the two groups, it means only that the number of differential genes compared to the mass of genes that are not differentially expressed is not enough to cluster the groups, which may be the case if the percentage of genes changed is small. There are some published evidences that recurrent spontaneous abortion, pre-term birth and preeclampsia share common mechanisms. The observation here seems to indicate that the uterus of women with PTB is identical to the one of term pregnancies. From the clinical data, the authors acted wisely in taking PTB without hypertension. But it is not clear from the data that the women that had PTB had a recurrent occurrence in this (not mentioned in Table S1).

5. In addition, the absence of DEG genes visible (Supplementary Figure 1B) is presented against FDR and not p-value which is unusual and I think a bot too stringent. Even random samples should give some significant genes (one out of 20) Using FDR leads to no significant genes, which is not a complete surprise. However, it does not mean that individual genes are not relevant for the difference between the two situations analyzed. Another interesting approach should be to check for enrichment of pathways using GSEA approaches, that do not rest on the establishment of thresholds that are always arbitrary. In sum excluding totally the existence of DEG between the two groups seem a bit rapid.

6. In the legend of Suppl figure S1, I do not understand 'Plot based on 728 genes labeled as fc (blue) and 17,748 genes labeled as none (purple)'. Nevertheless, it is admissible that the differences due to preterm or term are negligible compared to the differences induced by sPE.

The decomposition of the samples between sets for simulations is unusual in the field, but mathematically sound. The authors discovered 859 DEG at a threshold of 2-fold, and some genes (9) were validated by qRT-PCR.

7. It would be nice to analyze the data not only relative to a threshold of induction, but taking the complete dataset and using GSEA-like approaches to see whether they are consistent with the gene clusters found by threshold-volcano plots analyses.

Previous studies by the authors evaluate gene expression differences in decidualization either from human endometrial stromal cells (hESCs) from women with sPE versus normal. The authors identified 18 genes differentially expressed in vivo and in vitro in sPE compared to control.

8. The authors do question the observation of having only 18 common genes between in vivo and in vitro, out of 129 or 859 (in vitro – previous study- and in vivo -present study- respectively), and base them upon cell composition that is much more complex in vivo, which is reasonable. However, in the search for a signature, which is one of the justifications of the work, it could mean that the in vivo dataset could contain more 'relevant' genes compared to a in vitro model. The comparison with the results of Wang, which is much better (and focused better on a simple cell model) on the single-cell transcriptome may belong rather to a discussion rather than a results part of the paper. The validation by qRT-PCR is good, but the choice of the genes leads to a very high correlation (2D), that may be due to the use of only one induced gene (and only 5 genes used). About the actin as a control housekeeping gene, it is not sure that it is the best reporter gene in the uterine context. Generally, it is advised to normalize against the geometric mean of two to four different reporter genes.

A selection of 166 highly deregulated genes (>4 fold) is then selected, and were shown to be enough to separate efficiently the samples, which is not surprising, given that this corresponds to a semi-supervised analysis from genes found differential between the two gene sets.

9. The identification of Estrogen/progesterone receptor is not a surprise, when uterus function is concerned. I would suggest that the authors complete their analyses using network analyses such as provided by the combination of Stringdb and Cytoscape. This would help to visualize the pivotal position of ESR1 and PGR more clearly, or maybe to find other important hubs that were overlooked in their current study. Also using GSEA or other tools on transcription factor binding sites databases, showing the actual involvement (through measuring the enrichment, and calculation of FDRs) of the Estrogen Responsive Element and Progesterone Responsive Element would be an essential element that must be shown to prove a genuine enrichment in these cascades

Reviewer #3 (Recommendations for the authors):

There is much debate on whether the risk of PE is higher/lower in women with male vs. female fetuses. Additionally, the maternal response to pregnancy differs based on the sex of the fetus. The manuscript should therefore report on this variable and whether it impacts the results of the analysis.

Please clarify whether biopsies were done for clinical or research purposes only.

Why did the authors choose to use hg19 as their reference genome and not the better annotated hg38?

Methods (line 321) – What RIN value was defined as appropriate for making RNAseq libraries.

I cannot find anywhere where DD is defined. Can the authors please clearly define what they mean when they use DD?

Introduction – PE also significantly contributes to maternal mortality, not just infant.

Figure 3B – what is the axis label?

Figure 4 – I believe the authors mean sPE (not Spe).

Data availability: authors state that the RNAseq data is available for download in GEO.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Disrupted PGR-B and ESR1 signaling underlies defective decidualization linked to severe preeclampsia" for further consideration by eLife. Your revised article has been evaluated by Kathryn Cheah (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

Essential revisions:

The authors failed to properly respond to key reviewer concerns, which were considered to be necessary changes by the consensus of reviewer's discussion. Thus, the manuscript cannot be further considered until these points raised by the reviewers are appropriately addressed.

1. The editor and reviewers understood the method used for fold change thresholding. No additional explanation was required. The method is still wrong, as post-FDR fold-change thresholding is the problem that was raised and cannot be used in a rigorous analysis. Thus, the authors should implement one of the previously given options so that this major problem is corrected: (i) get rid of fold-change thresholding completely, or (ii) use a method that does it, while appropriately controlling for FDR (e.g. TREAT).

2. Similarly, using multiple corrections for GO and KEGG analysis is not a suggestion but a necessity for rigorous analysis, especially when performing so many tests. If term redundancy is a major concern for the authors in the GO database, they are advised to look at the GOSlim framework which was put together to specifically respond to this concern. In any case, enrichment results without FDR control are non-reproducible and not statistically supported.

3. The authors have not responded to the requirement to upload all code/scripts generated for the study. They refer to the software repository (i.e. edgeR), but not the actual code for analysis. The specific scripts that were written to run all analyses need to be made available as well on GitHub or as a supplement, as per eLife policy.

https://doi.org/10.7554/eLife.70753.sa1

Author response

Reviewer #1 (Recommendations for the authors):

1. There are major concerns about the statistics in the analysis of the dataset.

a. Fold-change thresholding is used in multiple places in the analysis (e.g. line 94, 112, 128, etc.). This is a big problem, as it is known to lead to poor FDR control in the case of statistical tests not designed to take fold change into account (like edgeR used here). Thus, the practice of fold-change filtering without FDR control is considered problematic by statisticians (see PMID: 19176553). Thus, either the authors need to update their results after removing this problematic criteria to properly control for FDR, or all analyses should be rerun with a method that takes fold-change into account for FDR control (e.g. the TREAT method PMID: 19176553).

We thank the reviewer for their feedback that allows us to explain in more detail the statistical analysis performed in the manuscript. First, we would like to highlight that all the analyses were performed with a proper false-discovery rate (FDR) cutoff as the basic criterion. Specifically, we used the exactTest function (edgeR package) that is recommended to make comparisons on datasets with a single-factor design, allowing us to achieve a global exploratory approach. This function uses the p-value adjustment method FDR, and a cutoff of 0.05 (FDR<0.05) was applied to identify significantly differentially expressed genes (DEGs). Once we obtained the DEGs, we selected those with a high fold change (FC≥2) to identify differences large enough to be biologically meaningful. We specified these two criteria (FDR<0.05 and FC≥2) in each line where the DEGs were described (e.g., Lines 94, 112, 128) and in the “Differential expression analysis” section in the Material and Methods (Lines 341–342), but we now provide more detail on the statistical analysis information used in this work.

– Specifically, we rewrote the “Quality control and pre-processing data” and “differential expression analysis” sections of the Materials and methods. These sections were replaced by a new section titled “RNA-seq analysis” (Lines 144-167): “Reads were mapped to the hg19 human genome transcriptome using the STAR (version 2.4.2a) read aligner (1). FastQC (version 0.11.2) was used to determine the quality of FASTQ files. The manipulation of SAM and BAM files was done with the software SAMtools (version 1.1) (2). To count the number of reads that could be assigned to each gene, we used HTSeq (version 0.6.1p1) (3) and BEDtools software (version 2.17.0) (4) to obtain gene coverage and work with bedFiles. Quality control filters at each program were used following the software package recommendations, and reads were filtered by mapping quality greater than 90%. Transcriptomic data were deposited in the Gene Expression Omnibus database (accession number GSE172381). The Bioconductor package edgeR (version 3.24.3) (5) was used to analyze differentially expressed genes. The trimmed mean of M-values normalization method was applied to our gene expression values. The exactTest function was used to find differentially expressed genes between groups. The p-value adjustment method was FDR with a cut-off of 0.05 (FDR<0.05). Once p-value was adjusted, significant deregulated genes with log2-fold-change ≥1 (FC≥2) were selected to perform gene ontology analysis and to formulate the signature encoding defective decidualization. edgeR analysis was carried out in R version 3.5.1. A volcano plot was created to visualize DEGs. For a better overview, we distinguished significant (FDR<0.05) and not significant (FDR≥0.05) DEGs with a high (FC≥2) or low fold-change (threshold FC<2).”

– We improved the explanation about the statistical analyses used in the Results section (Line 270-273): “Transcriptional analysis in the training set was performed by comparing gene expression patterns in sPE (n=17) and controls (n=12). This comparison revealed 859 significantly differentially expressed genes (DEGs) based on FDR<0.05 and at least two-fold differential expression between groups (FC≥2)” and (Line 293-295): “Eighteen genes were similarly differentially expressed with at least a two-fold change between groups (FDR<0.05 and FC≥2) (Figure 2A).”

– The legend of Figure 1B was also updated (Lines 673-675): “not significant-lowFC (FDR≥0.05, FC<2); not significant-highFC (FDR≥0.05; FC≥2); significant-lowFC (FDR<0.05; FC<2); and significant-highFC (FDR<0.05; FC≥2).”

– Supplementary file 3 has changed its name by “Figure1-Source data 1” and was included p-values, FDR values (corrected p-values), log FC, and FC, and we included how the criteria were applied in the legend, (lines 779-780) “Figure1-Source data 1. The 859 statistically differential expressed genes (FDR<0.05) with at least two-fold change (FC ≥ 2) in sPE vs control cases obtained from RNA-seq analysis.”

– The Supplementary file 4 have changed its name by “Figure 3-Source data 1” and legend was updated (783-784): “List of genes selected as the defective decidualization signature in sPE (166 DEGs with an FDR of <0.05 and an FC of ≥4).”

b. For the GO enrichment analysis (line 133), the authors only indicate a p-value filter (p < 0.005). Was there not implementation of multiple testing correction? The data should be re-analyzed after proper multiple hypothesis correction control (native DEseq2 FDR, Bonferroni, or other). In addition, what was the gene list background used to compute enrichment? The nature of the background has been shown to hugely impact results and needs to be explicitly detailed in the methods.

For this analysis, we followed the recommendation provided by the edgeR user’s guide:

“The p-values returned by goana and kegga are unadjusted for multiple testing. The authors have chosen not to correct automatically for multiple testing because GO terms and KEGG pathways are often overlapping, so standard methods of p-value adjustment may be very conservative. Users should be aware though that p-values are unadjusted, meaning that only very small p-values should be used for published results.”

For this reason, we selected a p-value filter of <0.005. Supplementary file 5 includes the p-value for each GO reported to provide statistical information to readers that are interested in using this information for their own research. Moreover, GO annotations emphasized in the text and in Figure 3 have a p-value of <0.005, with the vast majority lower than 0.00005. In addition, we agree with the reviewer that more information regarding the gene list background used to compute the analysis needs to be explicitly detailed, and we have added this information in Materials and methods.

– We rewrote the Enrichment analysis section of Materials and methods (Lines 180-188):

“GO analyses were conducted to obtain biological processes using the goana function in edgeR (62). The input genes were those 166 included in the fingerprinting (Figure 3-Source data 1) P-values returned automatically by the goana function are unadjusted for multiple testing because GO terms are often overlapping and standard methods of p-value adjustment may be very conservative. Thus, we kept GO terms with p-values of <0.005 (Figure3-Source Data 2); includes the p-value for each GO reported.”

– Supplementary file 5 have changed its name by “Figure 3-Source data 2”. The legend of was revised (785-787):

“Figure 3-Source data 2. Biological process GO terms computed by the 166 genes included in the DD fingerprint in sPE (N, number of genes associated with GO term; DE, number of genes differentially expressed in this GO term; P.DE, p-value).”

c. The authors reports using Student t-tests (section starting line 404). Unless normality is proven (e.g. with a Shapiro-Wilkes test), parametric tests that relie on normal distributions should never be used. We recommend that all tests be rerun with non-parametric equivalents (e.g. Wilcoxon or Mann Whitney) to guarantee these assumptions are not violated.

We followed the reviewer’s recommendation, and all the tests were rerun using the non-parametric method Wilcoxon to guarantee that assumptions of normality are not violated.

– The statistical analysis section of the Materials and methods was updated (Lines 232-234): “Clinical data were evaluated by Wilcoxon test for comparisons between sPE and control samples.”

– Supplementary file 1 was updated.

2. With human sample studies, it is unlikely that all samples were processed exactly in parallel (i.e. all samples collected the same day, kept in stabilizer for the same amount of time, etc.), which all but guarantees the existence of batch effects. There are also a number of stated parameters (i.e. age of the donor, time since pregnancy 1-8 years, day in cycle 22-32, etc.) that could influence results without a relevant biological driver.

a. In addition to an aggregated table (Table 1), the authors should generate a table for each anonymized sample for key characteristics so that the resource can be analyzed with these caveats considered both in this current paper and to facilitate all future reanalysis attempts. This should include for each sample individually:

i. age of the donor.

ii. time since last pregnancy.

iii. day of cycle for biopsy.

iv. date of processing/batch (including biopsy batch, library batch AND sequencing batch).

v. RIN of RNA sample.

vi. Sequencing depth (raw reads per library).

vii. Percent reads mapping to the genome reference (here hg19).

viii. Ethnicity of donor.

According to the reviewer’s suggestion, we have generated a Supplementary file 2 that includes the key characteristics indicated about each anonymized sample. The impact of these variables was considered in our analysis following the reviewer’s recommendations.

New Supplementary file 2 has been included in the manuscript (Line 112-114).

b. The authors should implement a multivariate model to account for all these potential confounding effects so as to focus on the core biological signal. In addition, the use of a software like SVA or RUVseq is recommended to remove unwanted technical variation prior to analyses.

We agree with the reviewer about the relevance of potential confounding factors of key characteristics on our transcriptomics results. We have applied an experimental design to minimize the impact of the main sources of batch effect on our gene expression data. Endometrial biopsies were preserved with stabilization solution and frozen at −80º until analysis to obtain the best-quality RNA. We used a balanced batch-group design. Specifically, case and control samples were included and were balanced and processed at the same time in the same technical batch, including the RNA extraction, library preparation, and sequencing. The batches performed were reduced to the minimal number of three. Additionally, before conducting the transcriptomic analysis, the confounding effects of main technical and clinical variables were tested. Once confounding effects were discarded, we performed a differential gene expression analysis. We included in the revised manuscript this additional information used as a basic starting point in our data analysis. In (Author response image 1), we provide to the reviewer some data from analyses to control for confounding effects in our experimental design. To exclude confounding effects in the gene expression analysis, we applied a principal variance component analysis (PVCA) to fit a mixed linear model using biological and technical variables as random effects to estimate and partition the total variability.

Author response image 1

The results demonstrated that the effect of biological and technical variables evaluated individually are less than 0.05, and their confounding effects are negligible. By contrast, the variable group (control and sPE) is the most plausible variable compared with the rest of the individual or double-interaction between variables (less than 0.1).– Materials and methods (Lines 128-129): “cDNA libraries from total RNA samples (n=40) were prepared using an Illumina TruSeq Stranded mRNA sample prep kit (Illumina, San Diego, CA) following a balanced batch-group design.”

– Results (Lines 247-248): “Biological and technical variables for each donor were considered to discard confounding effects on the transcriptomic profile (Supplementary file 2).”

– Discussion (Lines 437-439): “Having controlled for confounding effects of biological and technical variables, we consider that this misclustering is consistent with the nature of decidualization and its inherent variability.”

3. Additional methodological information is needed for long-term reproducibility of analyses.

a. For reproducibility of code and analyses, all analytical code for this study should be deposited in a repository such as Github or made available as a Supplementary file.

We also think that the reproducibility of our analysis is a basic scientific focus. This detailed information is specified in the revised manuscript, including that all sequencing data are available for download from the Gene Expression Omnibus (GEO; GSE172381). The software packages and their versions used in the analyses here are provided in the text, and the code needed to use the packages is publicly available on GitHub or Bioconductor (https://bioconductor.org/packages/release/bioc/html/edgeR.html).

See point #1a.

b. Please include all version numbers for all used software (e.g. R, etc.) packages and R packages (e.g. edgeR, goana, etc.), as well as all command parameters where relevant. The same should be applied to annotation databases (i.e. GO used here), and if a version number doesn't exist, date of access should be provided.

Thank you for this observation. We have included this information in the Materials and methods. See point #1a.

c. There is no mention of the softwares (nor versions) used to trim RNA-seq reads, to map to the hg19 genome or to aggregate counts to genes. All used software should be clearly named and their use described with all necessary options for reproducibility.

Thank you for this observation. We have included this information in the Materials and methods. See point #1a.

d. Line 84-88/Figure 1A, were the subtypes of "control" samples (pre vs. full-term) equally distributed in training vs. testing? Indeed, although the authors argue that the subtype made no difference on transcriptional profiles (i.e. Sup Figure 1A PCA), this reviews sees a general (if not clean separation) on PC1, that would probably be enhanced if sources of spurious variation (see point #2) were taken into account. In any case, the data provided cannot be used to support the statement "demonstrated an absence of clustering" (line 80), since the statistical tests used are not meant to prove the null hypothesis but to reject it when possible (i.e. confusion between type I and II error). Please amend the analysis to reflect this.

We agree with the reviewer that this preterm versus full-term analysis should be better explained. Our null hypothesis is that there is no differential gene expression across the two subtypes of control samples (preterm versus full-term). We performed a differential gene expression analysis with a proper FDR cutoff for multiple testing comparing preterm versus full-term, revealing that there was no statistical evidence of transcriptomic changes associated with the gestational week at delivery in our set of samples. Thus, gene expression patterns do not provide evidence for rejecting the null hypothesis of “no difference in gene expression between women with preterm or full-term labor who had never had preeclampsia”. This result is shown in the volcano plot in Figure supplement 1A In addition to this plot, we include a PCA because it is a technique for reducing the dimensionality of datasets, increasing interpretability but at the same time minimizing information loss. Thus, we consider it interesting for readers to keep the PCA figure in the manuscript to visualize variance among samples. However, we rewrote this part of the results to clarify that there were no significant differences between control groups. Furthermore, we would like to highlight that the two control subtypes were represented in both subsets of analyzed samples.

We rewrote this part of the Results to clarify this statement (Lines 248-262):

“Controls included women who had a preterm birth with no signs of infection (n=8) and women who gave birth at full term with normal obstetric outcomes (n=8). Transcriptomic profiles were compared by differential expression analysis, revealing no significant changes in the endometrial transcriptome between preterm and term controls [false-discovery rate (FDR) ≥0.05] (Figure 1—figure supplement 1A). Principal component analysis (PCA) supported that there was no underlying pattern of distribution depending on gestational age at delivery (Figure 1—figure supplement 1B). Once we ruled out bias on controls, we randomly split samples into two cohorts, a training set (70%) and a test set (30%) (Figure 1A).”

e. cDNA library cannot be validated for RIN (line 330), as only total RNA samples can be. The entire method section needs to be completely rechecked and rewritten for accuracy to avoid this kind of errors.

According to this recommendation, the Materials and methods section has been reviewed and rewritten when accuracy was required.

The Materials and methods section was updated (Lines 134-136): “cDNA libraries were quantified using an Agilent D1000 ScreenTape in a 4200 TapeStation system (Agilent Technologies Inc, Santa Clara, CA). Libraries were normalized to 10 nM and pooled in equal volumes”. See point #1 for more changes in the Materials and methods section.

f. All used QC filters should be stated explicitly in the methods (see lines 338-342).

We agree with this suggestion, and we have included the information in the manuscript.

Information included in the Materials and methods (lines 148-150): “Quality control filters in each program were used following the software package recommendations, and reads were filtered by mapping quality greater than 90%.”

4. Since RT-qpCR is known to be less sensitive and more prone to normalizing biases (due to the choice of control genes which may vary themselves) than RNA-seq (which performs unbiased normalization at the transcriptome-wide level), it is unclear why validation was performed on the same samples that were processed by RNA-seq (line 375-389). Unless the authors want to include RT-qPCR data on an INDEPENDENT cohort of patients, these results are circular and shouldn't be included in a revised manuscript. If RT-qPCR on a new cohort is included, make sure to include information on which normalizing amplicons were used for δ CT calculation.

We agree with the reviewer about the less sensitive value of RT-qPCR than that of RNA-seq. Following this reviewer’s suggestion, we removed the RT-qPCR data to validate our RNA-seq results from the same cohort of patients.

We removed these data in the revised manuscript. Specifically, data were removed from the Results (Lines 285-287; Lines 298-300). The Materials and methods section was updated (Lines 201-205) “Gene expression of IHH, MMP9, MSX2, ESR1, and PGR isoforms in the endometrial tissue from a subset of women with prior sPE (n=13) and controls (n=9) was obtained by RT-qPCR. Specific primers for each gene are described in Supplementary file 3.” Figure 2C, Figure 2D, and Figure supplement 2 were removed, and Supplementary file 3 was updated.

5. Can the authors explain why 3 specific samples (C20, 21, 22; line 147-152) do NOT cluster appropriately with their "signature"? Did they have any biological specificities (see need for sample-by-sample information as raised in point #2)? The fact that 3 control samples cluster with the sPE samples somewhat invalidate the signature as a signature of sPE, and may be the result of improper correction for batching or other technical (or irrelevant biological) noise. Please discuss explicitly what may be happening here in a revised manuscript.

We appreciate this comment. The new Supplementary file 2 (containing sample-by-sample information raised in point #2) shows that there was no biological or technical variable that allowed for the identification of the potential source for this misclustering. We suspect that this was due to the complexity of decidualization biology, which is explained in the Discussion section (lines 274–282). Decidualization is a highly dynamic process that shows interindividual variability. In addition, endometrial maturation starts around spiral arteries and extends to the entire endometrium. Consequently, random spatial sampling could affect the proportions of the decidualized cell types during biopsy collection. In fact, other authors, such as Munchel et al., 2020 (PMID: 32611681), have reported a misclustering of controls and PE cases using other approaches. This result is expected due to the high variability inherent to the pregnant human population, PE patients, and the endometrial cycle.

We decided to improve the explanation of this misclustering in the Discussion. Details related to Supplementary file 2 have been included in the Discussion section (Lines 437-439):

“Having controlled for confounding effects of biological and technical variables, we consider that this misclustering is consistent with the nature of decidualization and its inherent variability.”

Reviewer #2 (Recommendations for the authors):

In sum, I feel that this paper brings some interesting insights on the decidua status and preeclampsia; the degree of novelty is nevertheless unclear to me, and the clinical application seems far-fetched. There are several points that may be better presented and discussed.

Detail of the recommendations to the authors:

1. The study of the decidual transcriptome in preeclampsia has also already been performed by other teams in the past, rather at the moment of the disease than later. For instance, the paper of Mari Loset (Am J Obst Gynecol, 2011), is, I feel, an important base for comparison with the results presented here. So much so, that it used the same type of expressional approaches, and detected differential transcripts and cascades. I feel that the authors should analyze their own datasets in light of this type of seminal papers and discuss the commonalities and differences found.

Thank you for your comment. Loset’s work provided interesting insights of genetic canonical pathways and gene–gene interaction networks in decidua basalis from preeclamptic pregnancies at the end of pregnancy, as we previously did [Garrido-Gómez T. et al., 2017 (PMID: 28923940)]. The lack of coincidence between the transcriptomes of decidua during the second and third trimesters of pregnancy and decidualized endometrium in the secretory phase was expected. At the time of the decidua sampling in Loset’s paper, decidualization has already occurred, the decidua has been formed, and PE has been clinically manifested. Therefore, the pattern of gene expression in this tissue may be very different compared to cyclic decidualization without pregnancy or disease. However, following the reviewer’s suggestion, we compared our 859 DEGs with those 455 DEGs from Loset M et al., 2011. We found a modest overlap that included 28 genes, which is consistent with the expected outcome. We also compared our results with the dataset of decidua basalis and decidua parietalis dysregulated genes (79 and 227 DEGs, respectively) in sPE from our group [Garrido-Gómez T. et al., 2017 (PMID: 28923940)]. We found an overlap of 4 DEGs with decidua basalis and 14 DEGs with decidua parietalis. These results are consistent with Rabaglino MB et al., 2019, who performed multiple comparisons, including a comparison of the transcriptome of decidualized endometrial stromal cells from women who experienced sPE (PE-DEC) with the transcriptome of decidua basalis (PE-DB) and parietalis (PE-BP) from the placental bed of preeclamptic patients after delivery (PMID: 31356122). They found that there was little, if any, overlap of differentially expressed genes (DEGs) between PE-DEC and decidua specimens. Moreover, the decidual transcriptome at delivery differs from the decidual tissue of chorionic villi samples collected at ∼11.5 weeks. The authors explain, “These discordances may not be unexpected given the timing of the sample procurement—i.e., PE-DB (decidua basalis obtained from placental bed biopsy) and PE-BP (basal plate decidua obtained from delivered placentas) were procured during clinically active disease that by itself seems likely to perturb decidual gene expression, thus perhaps reflecting consequence rather than cause of disease.” Accordingly, our previous studies demonstrated impaired in vitro decidualization of ESC isolated from delivered placentas, consistent with the concept of endometrial antecedents of sPE (PMID: 28923940).

We agree with the authors that the lack of concordance between transcriptomes may be due to different type of samples (decidua and decidualized endometrium), and different timing of disease (samples at the end of the affected pregnancy versus samples years after) must show large differences in expressed genes. Thus, we did not include this type of comparison in the manuscript. Finally, in the manuscript, we include our own in vitro dataset because we wanted to test the abundance of new transcripts revealed for the first time in our in vivo bulk tissue approach. Also, we include the single-cell results because they were obtained via an in vivo approach using endometrial biopsies from healthy women across the endometrial cycle. Thus, these data could help identify which of the genes that were found to be dysregulated in women who experienced PE are potentially expressed by endometrial stromal cells. Altogether, these comparisons provide new data that we consider interesting due to its novelty.

Discussion (Lines 417-429):

“Previous reports of the decidual transcriptome in PE and sPE (10, 43), revealed the gene expression profile associated with the condition at the time of delivery. Here, we analyszed samples collected years after the affected pregnancy; thus, it could be interesting to find dysregulated genes in common among these approaches. We compared our transcriptomic results and the previously reported deysregulated decidual genes, obtaining 5.1% (10) and 6.2% (43) overlap in affected genes. These discordances may not be unexpected and are consistent with recent results (42). The decidua basalis transcriptome at delivery was compared with the transcriptome of decidua at ⁓11.5 gestational weeks and in vitro decidualized hESC from women who experienced sPE. Both analyses revealed little, if any, overlap between molecular signatures. In contrast, the signature encoding in vitro DD of hESC years after pregnancy overlapped significantly with the dysregulated profile found in decidual samples at the beginning of pregnancy (42). Thus, decidual gene expression patterns during the clinically active disease largely differ from those observed during endometrial decidualization at the end of the menstrual cycle and early pregnancy, perhaps reflecting a consequence rather than the origin of sPE.”

2. In the overrepresented pathways in the Loset's work, the regulation by ESR1 and PGR was not obvious. Does it mean that analyzing samples 3-4 years after the event (normal or pathological pregnancy) leads to remnants of the disease at the uterine level, or is there a genetic predisposition? This is a question that cannot be truly addressed by the present dataset, and is important if the objective is to define markers, as stated below.

The introduction gives the necessary details to understand the question raised, finding a signature that could help assessing the risk of preeclampsia.

Thank you for highlighting this interesting issue for discussion. In Loset’s work, the authors were focused on canonical pathways that could be affected by dysregulated genes in PE, which was valuable to increase the knowledge of PE pathogenesis. However, important differences between Loset’s work and our work need to be mentioned to answer this point. In Loset’s work, the authors did not examine in-depth the origin of the identified gene dysregulation. Biological processes or pathways tend to be quite specific, making it difficult to find obvious associations. We went one step further trying to investigate the association with hormone receptor activation due to the key role of ovarian hormones in decidualization. We tested this hypothesis, and this association became obvious from our results. However, Loset et al., did not address the same question.

Regarding the question “Does it mean that analyzing samples 3–4 years after the event (normal or pathological pregnancy) leads to remnants of the disease at the uterine level or is there a genetic predisposition?”, we agree with the reviewer that it could not be answered. Other analyses are needed, such as genetic variant analysis, to know if there is a genetic predisposition. We acknowledge the existence of this altered decidual transcriptome in sPE years after delivery, but there are other works that demonstrate that a decidualization defect occurs at both the end and the beginning of pregnancy (PMID: 28923940, PMID: 31356122). In addition, endometrial decidualization as a primary factor in pregnancy health and evidence to support the role of deficient decidualization in the origin of PE were presented last year (PMID: 32521725, PMID: 33007270). Therefore, genetic predisposition for a decidualization defect is plausible on this basis.

Discussion (Lines 498-501):

“However, the primary driver of the predisposition to undergo decidualization resistance and its link with the main risk factors of sPE remain to be determined. Our work is an important step toward the development of new strategies that enable early assessment of risk for sPE and might prompt new therapeutic strategies to treat this enigmatic pathological condition.”

3. However, the genetic part of preeclampsia is estimated by a heritability of ~55%. The genetic decomposition of this heritability presented by Cnattingius and coworkers (2004) indicate that only part of this genetics is connected to the maternal genetic background (which concerns two organs: the uterus and thus the decidua) and the placenta for half of its genetics. Part of the genetics of preeclampsia is associated only to the paternal genetics (estimated at ~1/3 of the heritability), therefore, it is clear that the risk cannot systematically be found by analysis of the female (uterine-decidua) expression profile only. Alternatively, the difference observed by the authors may be a consequence of the previous preeclampsia not based on a predisposition, but if this is the case, it limits considerably the usefulness of the markers found since preeclampsia, and in particular severe preeclampsia is at 75% a disease of the first pregnancy. The authors should recognize this as an important limit as an early prenatal screening strategy.

Preeclampsia is a syndrome in which different condition subtypes might coexist. Thus, to advance in its understanding, we have focused on the maternal contribution to severe preeclampsia through defective decidualization. Also, we are not claiming that our findings are related to heritability (DNA), but rather to the pathogenesis of the decidualization defect (RNA expression) involved in the origin of severe PE (PMID: 32521725; PMID: 33007270).

Discussion (lines 493-496):

“Preeclampsia is a syndrome in which different condition subtypes might coexist; here, we focused on the maternal contribution to sPE through DD. We are not claiming that our findings are related to sPE heritability (DNA), but instead, our findings are related to the pathogenesis of the decidualization defect (RNA expression) involved in the origin of PE (12, 22)”.

The initial analysis is based upon the endometrial tissue of 24 women that had a severe preeclampsia, 16 controls (8 preterm and 8 term births). Gestational age at delivery had apparently no influence on the 'control' groups. Since the samples were not collected at the end of gestation, which means that the occurrence of preterm may be connected to the placenta and not to the maternal uterine situation.

4. Nevertheless, there are only 8+8 samples so it is not possible to prove that in some case a decidual expression is not involved. In addition, using PCA to check that there is no clustering is not enough to certify that there are no differentially expressed genes that separate the two groups, it means only that the number of differential genes compared to the mass of genes that are not differentially expressed is not enough to cluster the groups, which may be the case if the percentage of genes changed is small. There are some published evidences that recurrent spontaneous abortion, pre-term birth and preeclampsia share common mechanisms. The observation here seems to indicate that the uterus of women with PTB is identical to the one of term pregnancies. From the clinical data, the authors acted wisely in taking PTB without hypertension. But it is not clear from the data that the women that had PTB had a recurrent occurrence in this (not mentioned in Table S1).

We thank the reviewer for the opportunity to improve the clarity of this point in our manuscript. We performed a gene expression analysis comparing preterm with full-term labor, revealing that there were no transcriptomic changes associated with the gestational age at delivery in our set of samples. This result is shown in the volcano plot in Figure supplement1A. In addition to this plot, we include a PCA because it is a widely used conventional method for reducing data complexity and visualizing the variance among samples. Thus, gene expression patterns do not provide evidence for rejecting the null hypothesis of “no difference in gene expression between women with preterm or full-term labor who had never had preeclampsia”. We agree with the reviewer that the sample size is limited to state conclusions such as “the uterus of women with PTB is identical to the one of term pregnancies”, but this is not our objective nor our intention. Our hypothesis is to demonstrate that patients that experience previous sPE have transcriptomic changes in the endometrium during the late secretory phase. To find the altered transcripts, we compared a case group with a control group, which was heterogeneous for the gestational age at delivery; but, testing differences due to preterm or full-term labor are negligible in controls compared to the differences induced in the cases that experience sPE.

We have rewritten this part of the Results to clarify this statement (Lines 248-256):

“Controls included women who had a preterm birth with no signs of infection (n=8) and women who gave birth at full term with normal obstetric outcomes (n=8). Transcriptomic profiles were compared by differential expression analysis, revealing no significant changes in the endometrial transcriptome between preterm and term controls [false-discovery rate (FDR) ≥0.05] (Figure 1—figure supplement 1A). Principal component analysis (PCA) supported that there was no underlying pattern of distribution depending on gestational age at delivery (Figure 1—figure supplement 1B).”

5. In addition, the absence of DEG genes visible (Supplementary Figure 1B) is presented against FDR and not p-value which is unusual and I think a bot too stringent. Even random samples should give some significant genes (one out of 20) Using FDR leads to no significant genes, which is not a complete surprise. However, it does not mean that individual genes are not relevant for the difference between the two situations analyzed. Another interesting approach should be to check for enrichment of pathways using GSEA approaches, that do not rest on the establishment of thresholds that are always arbitrary. In sum excluding totally the existence of DEG between the two groups seem a bit rapid.

The p-value adjustment method for multiple comparisons was an FDR with a cutoff of 0.05. Thus, we used the p-value adjusted for this analysis, and we applied the same method for the gene expression analysis between controls and sPE. We are comparing more than 18,000 genes between groups. With any p-value adjustment, we will obtain significant differences, but many of them would be false positives just by pure random chance. In this regard, FDR retains more significant p-values while increasing non-significant p-values than the Bonferroni method. Thus, we selected a less conservative adjustment method to identify differences if they exist, maintaining the FDR at 5%. In fact, the most preferable approach is controlling FDR as it not only reduces false positives but also minimizes false negatives (more details can be found in PMID: 30124010).

6. In the legend of Suppl figure S1, I do not understand 'Plot based on 728 genes labeled as fc (blue) and 17,748 genes labeled as none (purple)'. Nevertheless, it is admissible that the differences due to preterm or term are negligible compared to the differences induced by sPE.

The decomposition of the samples between sets for simulations is unusual in the field, but mathematically sound. The authors discovered 859 DEG at a threshold of 2-fold, and some genes (9) were validated by qRT-PCR.

Labels show the two criteria that we used to define the differentially expressed genes: adjusted p-value (FDR<0.05) and fold change (FC≥2). Blue dots represent genes that were not statistically different in expression, but that demonstrated an expression change (FDR≥0.05 and FC ≥2). Purple dots show those genes that were not significantly different in expression, with only a slight shift in expression (FDR≥0.05 and FC<2)

We rewrote the legend of Figure supplement 1 (752-758): “Figure 1—figure supplement 1. Transcriptomic analysis of control samples based on gestational age at delivery. (A) Volcano plot showing that there were no significant DEGs between controls according to gestational age at delivery. Labels show the two criteria that we used to define the differentially expressed genes: Adjusted p-value (FDR<0.05) and fold change (FC≥2). The plot is based on 728 genes labeled as “Not significant-FC≥2” (blue) and 17,748 genes labeled as “ Not significant-FC<2” (purple). (B) A PCA based on 18,476 genes, after removing genes with low expression, does not demonstrate clustering based on gestational age.”

7. It would be nice to analyze the data not only relative to a threshold of induction, but taking the complete dataset and using GSEA-like approaches to see whether they are consistent with the gene clusters found by threshold-volcano plots analyses.

We did not use an analysis of the data only relative to a threshold of induction, we took the complete dataset and used the edgeR package to obtain those genes that were significantly differentially expressed between groups (FDR<0.05). Once those DEGs were obtained, we emphasized the high fold change between groups (FC≥2 and FC≥4).

For the first filtering step, we used the exactTest function (edgeR package), which is recommended to make comparisons on datasets with a single-factor design, allowing us to achieve a global exploratory approach. This function uses the p-value adjustment method FDR (“False Discovery Rate”), and a cutoff of 0.05 was applied to identify DEGs. Once we obtained the significantly differentially expressed genes, we applied a second filtering step, fold change, to select those genes where differences were large enough to be biologically meaningful. To show visually the proportion of genes that pass through the filters (after analyzing the complete dataset), we used a volcano plot. We specified these two filters in each line, and we refer to DEGs and in the section “Differential expression analysis” from the Material and Methods, but we improve the information provided about the statistical analysis.

– Specifically, we rewrote “Quality control and pre-processing data” and “differential expression analysis” sections of the Materials and methods. These sections were replaced by a new section titled “RNA-seq analysis” (Lines144-167):

“Reads were mapped to the hg19 human genome transcriptome using the STAR (version 2.4.2a) read aligner (25). FastQC (version 0.11.2) was used to determine the quality of FASTQ files. The manipulation of SAM and BAM files was done with the software SAMtools (version 1.1) (26). To count the number of reads that could be assigned to each gene, we used HTSeq (version 0.6.1p1) (27) and BEDtools software (version 2.17.0) (28) to obtain gene coverage and work with bedFiles. Quality control filters in each program were used following the software package recommendations, and reads were filtered by mapping quality greater than 90%.. Transcriptomic data were deposited in the Gene Expression Omnibus database (accession number GSE172381). The Bioconductor package edgeR (version 3.24.3) (29) was used to analyze differentially expressed genes. The trimmed mean of M-values normalization method was applied to our gene expression values. The exactTest function was used to find differentially expressed genes between groups. The p-value adjustment method was FDR with a cut-off of 0.05 (FDR<0.05). Once p-value was adjusted, significant deregulated genes with log2-fold-change ≥1 (FC≥2) were selected to perform gene ontology analysis and to formulate the signature encoding DD. edgeR analysis was carried out in R version 3.5.1. A volcano plot was created to visualize DEGs. For a better overview, we distinguished significant (FDR<0.05) and not significant (FDR≥0.05) DEGs with a high (FC≥2) or low fold-change (threshold FC<2).”

– We improved the explanation about the statistical analysis applied in the Results section (Line 270-273):

“Transcriptional analysis in the training set was performed by comparing gene expression patterns in sPE (n=17) and controls (n=12). This comparison revealed 859 differentially expressed genes (DEGs) based on FDR<0.05 and at least two-fold differential expression between groups (FC≥2)” and (Line 293-295): “Eighteen genes were similarly differentially expressed with at least a two-fold change between groups (FDR<0.05 and FC ≥2) (Figure 2A).”

– The legend of Figure 1B was also updated (lines 672-674):

“not significant-lowFC (FDR≥ 0.05, FC<2); not significant-highFC (FDR≥0.05; FC≥2); significant-lowFC (FDR<0.05; FC< 2); significant-highFC (FDR<0.05; FC≥2).”

– Supplementary file 3 changed its name by Figure 1-Source data 1 and p-value, FDR values (corrected p-value), log FC, and FC, and we included how the criteria were applied in the legend (779-782):

“Figure 1-Source data 1. The 859 statistically differentially expressed genes (FDR<0.05) with at least two-fold change (FC≥2) in sPE versus control cases obtained from RNA-seq analysis.”

– The Supplementary file 4 changed its name by Figure 3- Source data 1 and legend was updated (lines 783-784):

“Figure 3- Source data 1. List of genes selected as the defective decidualization signature in sPE (166 DEGs with an FDR of <0.05 and an FC of ≥4).”

Previous studies by the authors evaluate gene expression differences in decidualization either from human endometrial stromal cells (hESCs) from women with sPE versus normal. The authors identified 18 genes differentially expressed in vivo and in vitro in sPE compared to control.

8. The authors do question the observation of having only 18 common genes between in vivo and in vitro, out of 129 or 859 (in vitro – previous study- and in vivo -present study- respectively), and base them upon cell composition that is much more complex in vivo, which is reasonable. However, in the search for a signature, which is one of the justifications of the work, it could mean that the in vivo dataset could contain more 'relevant' genes compared to a in vitro model. The comparison with the results of Wang, which is much better (and focused better on a simple cell model) on the single-cell transcriptome may belong rather to a discussion rather than a results part of the paper. The validation by qRT-PCR is good, but the choice of the genes leads to a very high correlation (2D), that may be due to the use of only one induced gene (and only 5 genes used). About the actin as a control housekeeping gene, it is not sure that it is the best reporter gene in the uterine context. Generally, it is advised to normalize against the geometric mean of two to four different reporter genes.

A selection of 166 highly deregulated genes (>4 fold) is then selected, and were shown to be enough to separate efficiently the samples, which is not surprising, given that this corresponds to a semi-supervised analysis from genes found differential between the two gene sets.

The comparison with the results of Wang (PMID: 32929266) was included in the Results section because this dataset was obtained from an in vivo approach using endometrial biopsies from healthy women across the endometrial cycle and identified gene expression patterns at a single-cell level. The analysis used in the current manuscript is from bulk tissue. Therefore, combining our dataset with the dataset from Wang provides insight into which of the genes that were identified as dysregulated in women who experienced preeclampsia may be expressed by endometrial stromal cells specifically. The data from this comparison analysis increases the specificity of our results.

We decided to remove the qRT-PCR data in the revised manuscript following the recommendation of Reviewer #1. The Reviewer suggested that RT-qPCR is less sensitive than RNA-seq and reported circular results because it was applied to the same cohort of patients.

The 166 highly deregulated genes were obtained comparing gene expression patterns between cases and controls from the training set. Then, these genes were applied to the test set of samples. This revealed an efficient separation of groups, validating the results observed in the training set.

The Materials and methods section was updated (Lines 201-205):

“Gene expression of IHH, MMP9, MSX2, ESR1, and PGR isoforms in the endometrial tissue from a subset of women with prior sPE (n=13) and controls (n=9) was obtained by RT-qPCR. Specific primers for each gene are described in Supplementary file 3.”

We removed these data in the revised manuscript. Specifically, data were removed from the Results (Lines 285-287; Lines 298-300). Figure 2C, Figure 2D, and Figure supplement 2 were removed, and Supplementary file 3 was updated.

9. The identification of Estrogen/progesterone receptor is not a surprise, when uterus function is concerned. I would suggest that the authors complete their analyses using network analyses such as provided by the combination of Stringdb and Cytoscape. This would help to visualize the pivotal position of ESR1 and PGR more clearly, or maybe to find other important hubs that were overlooked in their current study. Also using GSEA or other tools on transcription factor binding sites databases, showing the actual involvement (through measuring the enrichment, and calculation of FDRs) of the Estrogen Responsive Element and Progesterone Responsive Element would be an essential element that must be shown to prove a genuine enrichment in these cascades

Following the reviewer’s recommendation, the hub genes were determined by overlapping the top 20 genes obtained using two topological analysis methods of the cytoHubba plugin, maximal clique centrality (MCC) and maximum neighborhood component (MNC). In addition, genes responsive to estrogen and progesterone have been identified using the database of Human Transcription Factor Targets. Details are included in the revised manuscript.

– Materials and methods (Lines 195-198):

“Hub genes were extracted using the maximal clique centrality (MCC) and maximum neighborhood component (MNC) of the cytoHubba plugin (32). The overlapping genes identifiedy by the two topological analysis methods were selected as the hub genes.”

– Results (Lines 346-348):

“Regarding target genes of ER1 and PR, the database of Human Transcription Factor Targets (hTFtarget) reported 13 genes responsive to ER1 and 31 target genes modulated by PR, based on epigenomic, ChIP-seq, or motif evidence (36)”. Lines 358-361: “The hub genes were determined by overlapping the top 20 genes obtained using two topological analysis methods in the cytoHubba plugin (31), maximal clique centrality (MCC) and maximum neighborhood component (MNC). Sixteen genes were selected, all of which were downregulated.”

Reviewer #3 (Recommendations for the authors):

There is much debate on whether the risk of PE is higher/lower in women with male vs. female fetuses. Additionally, the maternal response to pregnancy differs based on the sex of the fetus. The manuscript should therefore report on this variable and whether it impacts the results of the analysis.

We included 24 women who experienced sPE, of whom 9 had a female fetus and 11 had a male fetus (data for two patients were unavailable). The PCA and heat map in (Author response image 2) represent this variable, showing a heterogeneous distribution of samples [female (yellow) and male (green) fetuses]. Thus, we do not include this result in our manuscript since we do not focus on sPE risk factors, and we the sex of the fetus in our limited series does not impact the results of our transcriptomic data.

Author response image 2
Validation of the DD fingerprint in sPE.

(A) PCA based on 166 genes included in the fingerprinting in the training set. Each sample is represented as a colored point (blue, control; orange, sPE). (B) Heatmap dendogram of expression of the 166 genes included in the final fingerprinting for each sample of the training set (control, n=12; sPE, n=17). Sex of the fetus is represented by color (yellow, female; green, male). (C) PCA based on the fingerprinting in the test set. Each sample is represented as a colored point (blue, control; orange, sPE). (D) Heatmap dendogram of expression of the 166 genes included in the final fingerprinting for each sample of the test set (control, n=4; sPE, n=7). Sex of the fetus is represented by color (yellow, female; green, male).

Please clarify whether biopsies were done for clinical or research purposes only.

Biopsies were performed for research purposes only. This study was approved by the Clinical Research Ethics Committee of Hospital La Fe (Valencia, Spain) (2011/0383).

We clarified this observation in the manuscript (Line 92-94):

“Endometrial samples were obtained for research purposes during late secretory phase in 24 women who had developed sPE in a previous pregnancy and in 16 women with no history of sPE with full term (n=8) and preterm pregnancies (n=8) as controls.”

Why did the authors choose to use hg19 as their reference genome and not the better annotated hg38?

The latest build of the human reference genome, commonly nicknamed hg38, greatly expanded the repertoire of ALT contigs. These represent alternate haplotypes and have a significant impact on our power to detect and analyze genomic variation that is specific to populations that carry alternate haplotypes. Thus, this version is strongly recommended for analyzing genetic variation. However, we were interested in analyzing gene expression patterns instead of sequence variation. The existence of different versions causes confusion, and part of the problem is that many bioinformatic tools fail to enforce consistent use of a specific reference. This allows the unwary user to switch reference genomes halfway through a project without realizing that their comparisons suddenly become worthless (because, for example, now all the positions are shifted by some coordinate index). Hg19 is still the most widely used version; consequently, many studies, and even clinical databases and bioinformatic tools, are based on it. Thus, version hg19 is worthwhile for gene expression analysis and facilitates the use of our dataset by the scientific community.

We included these details in the Materials and methods (lines 144-145):

“Reads were mapped to the hg19 human genome transcriptome using the STAR (version 2.4.2a) read aligner (25).”

Methods (line 321) – What RIN value was defined as appropriate for making RNAseq libraries.

RIN values ranged from 4.9 to 9.2. This detail has been included (Line 124).

I cannot find anywhere where DD is defined. Can the authors please clearly define what they mean when they use DD?

Thank you for this observation to clarify this abbreviation in the manuscript. We defined the term in Lines 63-65. The abbreviation is now included:

“Defective decidualization (DD) entails the inability of the endometrial compartment to undertake tissue differentiation, leading to aberrations in placentation and compromising pregnancy health (12)”.

Introduction – PE also significantly contributes to maternal mortality, not just infant.

We have included the PE contribution to maternal mortality in this sentence (line 34).

Change in line 49-51:

“PE is characterized by the onset of hypertension, proteinuria, and other signs of maternal vascular damage that contribute to maternal and neonatal mortality and morbidity (1)”.

Figure 3B – what is the axis label?

The axis label refers to “Enrichment index”, which is calculated by −log(p-value). In Figure 3B, the axis label has been included in the figure and is explained in the legend.

Figure 4 – I believe the authors mean sPE (not Spe).

Yes, we mean sPE. Thank you for the observation. Spe was replaced by sPE in Line 721

Data availability: authors state that the RNAseq data is available for download in GEO.

Transcriptomic data were deposited in the Gene Expression Omnibus database (accession number GSE172381), as specified in the Materials and methods (RNA-seq analysis). However, after this recommendation, we decided to include the accession number in the Results section.

GSE code was included in the Results section (Lines 238-240):

“To identify transcriptomic alterations during decidualization in sPE, we applied global RNA sequencing (RNA-seq) to endometrial biopsies obtained in the late secretory phase from women who developed sPE in a previous pregnancy (n=24) and from controls who never had sPE (n=16) (GSE172381).”

[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]

Essential Revisions (for the authors):

The authors failed to to properly respond to key reviewer concerns, which were considered to be necessary changes by the consensus of reviewer's discussion. Thus, the manuscript cannot be further considered until these points raised by the reviewers are appropriately addressed.

We would like to thank the reviewers and editors for their time and effort devoted to improving our manuscript. We addressed the key reviewer concerns in the revised manuscript, please find bellow our point-by-point answer.

1. The editor and reviewers understood the method used for fold change thresholding. No additional explanation was required. The method is still wrong, as post-FDR fold-change thresholding is the problem that was raised and cannot be used in a rigorous analysis. Thus, the authors should implement one of the previously given options so that this major problem is corrected: (i) get rid of fold-change thresholding completely, or (ii) use a method that does it, while appropriately controlling for FDR (e.g. TREAT).

Following the reviewer suggestion we used the method TREAT to test formally the hypothesis (with associated p-values) that the differential expression is greater than a given threshold. Changes have been implemented along the manuscript since the analysis was rerun with a different method (TREAT instead of ExacTest). Thus, we updated Materials and Methods, results and figures.

2. Similarly, using multiple corrections for GO and KEGG analysis is not a suggestion but a necessity for rigorous analysis, especially when performing so many tests. If term redundancy is a major concern for the authors in the GO database, they are advised to look at the GOSlim framework which was put together to specifically respond to this concern. In any case, enrichment results without FDR control are non-reproducible and not statistically supported.

According to this recommendation, multiple corrections for GO analysis was applied. We rewrote the part of GO analysis in Materials and Methods (lines 170-176) and results (lines 277-292).

3. The authors have not responded to the requirement to upload all code/scripts generated for the study. They refer to the software repository (i.e. edgeR), but not the actual code for analysis. The specific scripts that were written to run all analyses need to be made available as well on GitHub or as a supplement, as per eLife policy.

We made available the scripts that were written to run all analyses. They could be found in Github (https://github.com/mclemente-igenomix/garrido_et_al_2021). We included the sentence “Custom scripts are available on GitHub at https://github.com/mclemente-igenomix/garrido_et_al_2021” in lines 159, 166, and 176. Also, the link was included in the key resources table.

https://doi.org/10.7554/eLife.70753.sa2

Article and author information

Author details

  1. Tamara Garrido-Gomez

    Igenomix Foundation, INCLIVA, Valencia, Spain
    Contribution
    Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Validation, Visualization, Writing – original draft, Writing – review and editing
    Contributed equally with
    Nerea Castillo-Marco and Mónica Clemente-Ciscar
    For correspondence
    tamara.garrido@igenomixfoundation.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6584-4832
  2. Nerea Castillo-Marco

    Igenomix Foundation, INCLIVA, Valencia, Spain
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Validation, Writing – original draft
    Contributed equally with
    Tamara Garrido-Gomez and Mónica Clemente-Ciscar
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4817-4777
  3. Mónica Clemente-Ciscar

    Igenomix, Valencia, Spain
    Contribution
    Data curation, Formal analysis, Investigation, Writing – original draft
    Contributed equally with
    Tamara Garrido-Gomez and Nerea Castillo-Marco
    Competing interests
    No competing interests declared
  4. Teresa Cordero

    Igenomix Foundation, INCLIVA, Valencia, Spain
    Contribution
    Investigation, Methodology, Project administration, Validation
    Competing interests
    No competing interests declared
  5. Irene Muñoz-Blat

    Igenomix Foundation, INCLIVA, Valencia, Spain
    Contribution
    Investigation, Methodology, Writing – original draft
    Competing interests
    No competing interests declared
  6. Alicia Amadoz

    Igenomix, Valencia, Spain
    Contribution
    Data curation, Formal analysis, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3915-0404
  7. Jorge Jimenez-Almazan

    Igenomix, Valencia, Spain
    Contribution
    Data curation, Formal analysis, Investigation
    Competing interests
    No competing interests declared
  8. Rogelio Monfort-Ortiz

    Department of Obstetrics and Gynecology, University and Polytechnic La Fe Hospital, Valencia, Spain
    Contribution
    Investigation, Project administration, Resources
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7931-8609
  9. Reyes Climent

    Department of Obstetrics and Gynecology, University and Polytechnic La Fe Hospital, Valencia, Spain
    Contribution
    Investigation, Project administration, Resources
    Competing interests
    No competing interests declared
  10. Alfredo Perales-Marin

    1. Department of Obstetrics and Gynecology, University and Polytechnic La Fe Hospital, Valencia, Spain
    2. Department of Obstetrics and Gynecology, School of Medicine, Valencia University, Valencia, Spain
    Contribution
    Resources, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2221-2560
  11. Carlos Simon

    1. Igenomix Foundation, INCLIVA, Valencia, Spain
    2. Department of Obstetrics and Gynecology, School of Medicine, Valencia University, Valencia, Spain
    3. Obstetrics & Gynecology, BIDMC Harvard University, Boston, United States
    Contribution
    Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Writing – review and editing
    For correspondence
    carlos.simon@uv.es
    Competing interests
    No competing interests declared

Funding

Carlos III Health Institute (Grant PI19/01659 (MCIU/AEI/FEDER, UE))

  • Tamara Garrido-Gomez

Spanish Generalitat Valenciana (PhD Student grant FDGENT/2019/008)

  • Nerea Castillo-Marco

Spanish Ministry of Science and Innovation (PhD Student grant PRE2019-090770)

  • Irene Muñoz-Blat

Spanish Ministry of Science and Innovation (Grant RTI2018-094946-B-100 (MCIU/AEI/FEDER, UE))

  • Carlos Simon

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are grateful to Dr. Alfredo Perales and the Obstetric Area from the Hospital University and Polytechnic La Fe Hospital for invaluable help in the enrollment of participants and obtaining the tissue samples that made this study possible. We thank the University and Polytechnic La Fe Hospital recruiters Rogelio Monfort, Reyes Climent, Laura Rubert, Joana Dasí, and Julia Escrig for their assistance in compiling clinical data. We are indebted to the patient participants. This work was supported by the grant PI19/01659 (MCIU/AEI/FEDER, UE) from the Spanish Carlos III Institute awarded to TG-G. NC-M was supported by the PhD program FDGENT/2019/008 from the Spanish Generalitat Valenciana. IM-B was supported by the PhD program PRE2019-090770 and funding was provided by the grant RTI2018-094946-B-100 (MCIU/AEI/FEDER, UE) from the Spanish Ministry of Science and Innovation with CS as principal investigator. This work was funded partially by Igenomix S.L.

Ethics

Human subjects: This study was approved by the Clinical Research Ethics Committee of Hospital La Fe (Valencia, Spain) (2011/0383), and written informed consent was obtained from all participants before tissue collection and all samples were anonymized (Included in Methods section - Human donors).

Senior Editor

  1. Kathryn Song Eng Cheah, The University of Hong Kong, Hong Kong

Reviewing Editor

  1. Bérénice A Benayoun, University of Southern California, United States

Version history

  1. Received: May 27, 2021
  2. Preprint posted: July 24, 2021 (view preprint)
  3. Accepted: September 23, 2021
  4. Version of Record published: October 28, 2021 (version 1)

Copyright

© 2021, Garrido-Gomez et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 945
    Page views
  • 145
    Downloads
  • 8
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Tamara Garrido-Gomez
  2. Nerea Castillo-Marco
  3. Mónica Clemente-Ciscar
  4. Teresa Cordero
  5. Irene Muñoz-Blat
  6. Alicia Amadoz
  7. Jorge Jimenez-Almazan
  8. Rogelio Monfort-Ortiz
  9. Reyes Climent
  10. Alfredo Perales-Marin
  11. Carlos Simon
(2021)
Disrupted PGR-B and ESR1 signaling underlies defective decidualization linked to severe preeclampsia
eLife 10:e70753.
https://doi.org/10.7554/eLife.70753

Further reading

    1. Cell Biology
    2. Chromosomes and Gene Expression
    Maikel Castellano-Pozo, Georgios Sioutas ... Enrique Martinez-Perez
    Short Report Updated

    The cohesin complex plays essential roles in chromosome segregation, 3D genome organisation, and DNA damage repair through its ability to modify DNA topology. In higher eukaryotes, meiotic chromosome function, and therefore fertility, requires cohesin complexes containing meiosis-specific kleisin subunits: REC8 and RAD21L in mammals and REC-8 and COH-3/4 in Caenorhabditis elegans. How these complexes perform the multiple functions of cohesin during meiosis and whether this involves different modes of DNA binding or dynamic association with chromosomes is poorly understood. Combining time-resolved methods of protein removal with live imaging and exploiting the temporospatial organisation of the C. elegans germline, we show that REC-8 complexes provide sister chromatid cohesion (SCC) and DNA repair, while COH-3/4 complexes control higher-order chromosome structure. High-abundance COH-3/4 complexes associate dynamically with individual chromatids in a manner dependent on cohesin loading (SCC-2) and removal (WAPL-1) factors. In contrast, low-abundance REC-8 complexes associate stably with chromosomes, tethering sister chromatids from S-phase until the meiotic divisions. Our results reveal that kleisin identity determines the function of meiotic cohesin by controlling the mode and regulation of cohesin–DNA association, and are consistent with a model in which SCC and DNA looping are performed by variant cohesin complexes that coexist on chromosomes.

    1. Chromosomes and Gene Expression
    2. Developmental Biology
    Airat Ibragimov, Xin Yang Bing ... Paul Schedl
    Research Article Updated

    Though long non-coding RNAs (lncRNAs) represent a substantial fraction of the Pol II transcripts in multicellular animals, only a few have known functions. Here we report that the blocking activity of the Bithorax complex (BX-C) Fub-1 boundary is segmentally regulated by its own lncRNA. The Fub-1 boundary is located between the Ultrabithorax (Ubx) gene and the bxd/pbx regulatory domain, which is responsible for regulating Ubx expression in parasegment PS6/segment A1. Fub-1 consists of two hypersensitive sites, HS1 and HS2. HS1 is an insulator while HS2 functions primarily as an lncRNA promoter. To activate Ubx expression in PS6/A1, enhancers in the bxd/pbx domain must be able to bypass Fub-1 blocking activity. We show that the expression of the Fub-1 lncRNAs in PS6/A1 from the HS2 promoter inactivates Fub-1 insulating activity. Inactivation is due to read-through as the HS2 promoter must be directed toward HS1 to disrupt blocking.