Transcriptomics-Based Subphenotyping of the Human Placenta Enabled by Weighted Correlation Network Analysis in Early-Onset Preeclampsia With and Without Fetal Growth Restriction

Background: Placental disorders contribute to pregnancy complications, including preeclampsia and fetal growth restriction (FGR), but debate regarding their specific pathobiology persists. Our objective was to apply transcriptomics with weighted gene correlation network analysis to further clarify the placental dysfunction in these conditions. Methods: We performed RNA sequencing with weighted gene correlation network analysis using human placental samples (n=30), separated into villous tissue and decidua basalis, and clinically grouped as follows: (1) early-onset preeclampsia (EOPE)+FGR (n=7); (2) normotensive, nonanomalous preterm FGR (n=5); (2) EOPE without FGR (n=8); (4) spontaneous idiopathic preterm birth (n=5) matched for gestational age; and (5) uncomplicated term births (n=5). Our data was compared with RNA sequencing data sets from public databases (GSE114691, GSE148241, and PRJEB30656; n=130 samples). Results: We identified 14 correlated gene modules in our specimens, of which most were significantly correlated with birthweight and maternal blood pressure. Of the 3 network modules consistently predictive of EOPE±FGR across data sets, we prioritized a coexpression gene group enriched for hypoxia-response and metabolic pathways for further investigation. Cluster analysis based on transcripts from this module and the glycolysis/gluconeogenesis metabolic pathway consistently distinguished a subset of EOPE±FGR samples with an expression signature suggesting modified tissue bioenergetics. We demonstrated that the expression ratios of LDHA/LDHB and PDK1/GOT1 could be used as surrogate indices for the larger panels of genes in identifying this subgroup. Conclusions: We provide novel evidence for a molecular subphenotype consistent with a glycolytic metabolic shift that occurs more frequently but not universally in placental specimens of EOPE±FGR.

P lacental function within physiological limits is essential for successful pregnancy outcomes. Abnormal development or dysfunction of this vital organ can lead to severe pregnancy complications, inclusively referred to as pregnancy-related placental syndromes. 1,2 Disorders associated with disordered placentation include fetal growth restriction (FGR) and preeclampsia, 2 both of which increase the acute risk of perinatal morbidity and mortality. 1,3,4 These disorders can predispose pregnant individuals and their offspring to chronic health problems. 1 Preeclampsia has traditionally been defined clinically by new-onset or worsening hypertension during pregnancy >20 weeks with accompanying proteinuria or other signs of organ system involvement. 5 FGR, which can accompany preeclampsia or occur independently, is typically detected prior to birth through ultrasonographic estimation of fetal weight. When it occurs, placental insufficiency is strongly implicated in the pathophysiology of FGR and preeclampsia alike; however, the precise etiological underpinnings of these complex conditions remain incompletely understood. 2,3,6 Final common pathway models for both disease states have been proposed. These posit that diverse causes can converge to produce suboptimal uteroplacental perfusion, which in turn may manifest itself in stereotyped, albeit inconsistent, clinical findings. 6 Weighted gene correlation network analysis (WGCNA) is a data reduction technique applied widely in gene expression studies to identify and functionally categorize clusters of highly correlated transcripts (coexpression modules). 7 The resulting clusters can then be related to clinical variables (eg, blood pressure) using consolidation metrics such as eigengenes 8 and metagenes. 9 Since modules thus constructed are impartial to clinical outcome, WGCNA can be applied when analyzing data sets with unknown sample stratifications. Translational examples include use in identifying tumor subgroups with differing outcomes but common presentations. 10 By distinguishing molecular subtypes and their links to candidate prognostic markers and druggable targets, this approach allows investigators to crystalize transcriptome expression data in a potentially actionable way. 7 Recognizing the considerable heterogeneity in outcomes for patients presenting with preeclampsia and FGR, Leavey et al 11 and Gibbs et al 12 used unsupervised clustering approaches for class discovery using gene expression profiling in a large cohort of placental samples from affected pregnancies. Their data-driven classifications revealed several robust subclusters that differed in clinical and histopathologic findings as well as molecular-level functional enrichments. These studies illustrated the potential for functional genomics in molecular phenotyping as an adjunct to clinical diagnostic criteria in classifying preeclampsia and FGR and offered new insights into the pathoetiology of these conditions. Inspired by these findings, we used WGCNA in combination with transcriptome profiling by RNA sequencing (RNA-seq) to (1) define other possible placental subphenotypes in preeclampsia and FGR; and (2) further elucidate dysregulated placental processes associated with these syndromes.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request. The RNA-seq data generated for this publication have been deposited in NCBI's Gene Expression Omnibus (GSE203507). Detailed Methods are provided in the Supplemental Material.

Study Approval
This study was approved by the Institutional Review Board of Yale University, The Ohio State University, the Abigail Wexner

Novelty and Significance
What Is New?
This study demonstrates the following: A subphenotype consistent with a glycolytic metabolic shift occurs more frequently but not universally in placentas from pregnancies complicated by preeclampsia. The ratios of LDHA/LDHB and PDK1/GOT1 mRNA transcripts may serve as surrogate indices for larger gene panels to identify metabolically stressed placentas.

What Is Relevant?
The 2 placental transcript ratios described in this study may add additional precision to the molecular subphenotyping of the placenta beyond want could be inferred from clinical characteristics alone.

Clinical/Pathophysiological Implications?
Our findings add to general knowledge that preeclampsia is a heterogeneous hypertensive syndrome, and this heterogeneity extends to the placenta.

RNA Extraction and Bulk RNA-seq
Villous tissue (VT) and basal plate decidua basalis (DB) were collected, and total RNA was isolated, as previously reported. 13,16 RNA-seq was performed using the Illumina HiSeq 2500 system.

Differential Abundance Analysis
Differential transcript abundance was determined using DESeq2 v1.32.0. Surrogate variable analysis was performed to control for unrecognized heterogeneity. Statistical models included clinical outcome, assigned fetal sex at delivery, and the first 2 components of the surrogate variable analysis as covariates. To minimize batch effects, the RNA-seq data were generated in 2 rounds of sequencing. The PTB and TB specimens reported previously [13][14][15] were isolated and sequenced concurrently with samples S16, S20, S21, and S22 (Table S2).

Weighted Gene Correlation Network Analysis
WGCNA was performed using the top quartile (≈4000) of normalized transcript counts exhibiting the greatest variability in each placental region (VT and DB). Correlated gene modules were identified from coexpression networks using the WGCNA R library.

Module-Trait Associations
Eigengene-based module-trait associations for continuous clinical variables were computed using the WGCNA R package. For categorical variables, we adopted a metagene-centered approach. Briefly, associations between binary outcomes and metagenes, defined as the arithmetic mean of normalized, variance-stabilizing transformed counts for selected transcripts, 9 were assessed using receiver operating characteristic (ROC) analysis. For each module, the metagene comprising the top 10 transcripts most positively correlated with select traits was divided by the metagene composed of the top 10 transcripts most negatively correlated with those same traits, yielding a metagene ratio (MGR).

Comparison with Previously Published Data Sets
To assess the consistency of our present findings with those of prior studies, we identified 3 placenta-derived RNA-seq data sets for further analysis: (1) GSE114691, 17 comprising 79 samples; (2) GSE148241, 18 consisting of 41 samples; and (3) PRJEB30656, 19 consisting of 10 samples (Table S3). For joint analyses, variance-stabilizing transformed counts were adjusted for study-dependent batch effects.

Functional Enrichment Analysis
We performed functional enrichment (hypergeometric testing) on prioritized groups of transcripts using clusterProfiler in R.

qPCR Cross-Validation
We selected 4 transcripts for quantitative polymerase chain reaction (qPCR) cross-validation using the general methodology described previously 16 : lactate dehydrogenase A (LDHA), lactate dehydrogenase B (LDHB), pyruvate dehydrogenase kinase 1 (PDK1), and glutamic-oxaloacetic transaminase 1 (GOT1). These transcripts were chosen based on their relevance to core placental metabolic pathways. [20][21][22][23] TaqMan gene expression assays (Thermo Fisher Scientific) were used for quantification (see the Major Resources Table). Relative abundance was calculated using the comparative Ct method. The geometric mean of the Ct values for β2-microglobulin (B2M) and ribosomal protein L30 (RPL30) was used as a reference in each reaction, performed in duplicate. Correlations between the expression levels from RNA-seq and qPCR experiments were examined using Pearson product-moment correlation (on log-transformed data).

Statistics
All statistical analyses were performed using a combination of functions and packages within the R v4.1.0, in addition to Prism v9 (GraphPad Software, La Jolla, CA) and MedCalc v20.027 (MedCalc Software Ltd, Ostend, Belgium). Clinical characteristics among the study groups were evaluated using ANOVA or Kruskal-Wallis testing (for continuous variables) and χ 2 tests (for categorical variables) and for these comparisons, a P <0.05 was considered statistically significant. Adjustments for multiple comparisons were performed using the Benjamini-Hochberg procedure with false discovery rate (FDR) <0.1 considered significant. Additional details are provided in the figure legends and Supplemental Material.

Overview of Study Population and Differential Abundance Analysis
We performed RNA-seq profiling of paired VT and DB placental samples grouped by clinical phenotype: (1) EOPE+FGR; (2) normotensive, nonanomalous preterm FGR; (3) EOPE without FGR; (4) spontaneous idiopathic PTB without FGR, matched for GA; and (5) uncomplicated TB. The 2 control groups (PTB and TB) were chosen to account for the changes in gene expression associated with spontaneous preterm delivery or advancing GA. The demographic, clinical, and pregnancy outcome data associated with these placental samples are summarized in Tables S1 and S2.
Similarity analysis with hierarchical clustering revealed that PTB and TB placental samples were the most similar groups overall, with EOPE and EOPE+FGR samples forming a more distantly related group and FGR occupying an intermediate position ( Figure 1). Consistently, by principal component (PC) analysis, we noted that PTB and TB placental samples were more closely aggregated than the EOPE, EOPE+FGR, and FGR samples ( Figure 1B and 1D). More detailed pairwise comparisons between the groups recapitulated these general observations ( Figure S1 through S3).

Correlation Network Analysis Reveals Placental Transcript Modules That Correlate With Clinical Variables Semi-Independently
We next generated weighted correlation networks to further assess potentially informative gene expression patterns in relation to select clinical variables. This approach enables unsupervised, data-driven assessment of the interplay between molecular signatures and various clinical variables. 7 Eight coexpression modules were identified in the VT samples (Figure 2A and 2B; Table S4). Moduletrait correlation analysis showed that, among VT samples, modules VT 2 (blue, 1001 transcripts) and VT 3 (black, 391 transcripts) were the most strongly correlated with blood pressure, whereas module VT 1 (orange, 1066 transcripts) showed the strongest correlation with birthweight ( Figure 2C). Modules VT 4 (turquoise, 327 transcripts) and VT 7 (brick, 238 transcripts) also exhibited a significant correlation with birthweight in addition to blood pressure, although VT 6 (gold, 269 transcripts) was correlated with birthweight only ( Figure 2C).
Functional enrichment analysis of the VT orange module demonstrated overrepresentation for proteasome pathway mRNAs (including transcripts encoding multiple proteasome subunits) and ribosomal protein subunit transcripts (all FDR <0.006, hypergeometric test; Figure 2D). The VT blue module included numerous preeclampsia -associated transcripts curated in the Comparative Toxicogenomics Database, such as ENG, FLT1, and LEP (Table S4), and was enriched (FDR <0.006) for several membrane receptor signaling pathways (Figure 2D). The VT black module was chiefly characterized (FDR <0.081) by transcripts associated with the HIF (hypoxia-inducible factor) pathway and glucose metabolism ( Figure 2D). The VT turquoise module showed overrepresentation for transcripts involved in extracellular matrix interactions and related signaling pathways (FDR <0.04; Figure 2D).
In DB tissues, 6 subnetworks of correlated expression were recognized ( Figure 3, Table S5). The DB mustard, DB tomato, and DB magenta modules were most highly correlated with birthweight, although the DB brown module was the only DB subnetwork with a significant blood pressure correlation ( Figure 3C). The DB mustard module ( Figure 3D) was enriched for expressed genes involved in growth factor signaling and nuclear envelope function (FDR <0.07, hypergeometric test). Transcripts of the DB tomato module ( Figure 3D) featured proteincoding transcripts involved in histone methylation, processing of mRNA, and Notch signaling (FDR <0.06). The DB magenta group was highly enriched (FDR <0.001) for genes involved in host immune response pathways ( Figure 3D). The DB brown module included transcripts encoding proteins involved in the HIF and transforming growth factor-β pathways in addition to those involved in the renin-angiotensin-aldosterone response, but these enrichments were not statistically significant (FDR >0.1; Figure 3D). The full module enrichment analysis results are shown in Figures S4 and S5.
Modules within each anatomical placental region shared no common transcripts by algorithmic design. Between these regions, the degree of similarity among individual modules measured using the Jaccard index ranged from moderate (Jaccard index, 17.2% between VT orange and DB tomato) to trivial (Jaccard index ≤3% for 30 module pairs; Figure S6).

Coexpression Network Module Reproducibility and Predictive Performance Across Studies
Between the 2 structural placental regions (VT versus DB), all transcript modules except for the DB mustard group showed statistical evidence for preservation (Figure 4A and 4B). Given this, we proceeded to examine network module reproducibility more generally, using 3 previously published placental RNA-seq data sets in which similar clinical phenotypes were studied: (1) GSE148241, 18 consisting of EOPE (n=3), EOPE+FGR (n=6), and term control (n=32) specimens; (2) GSE114691, 17 comprising 79 samples including FGR (n=18), EOPE (n=20), EOPE+FGR (n=20), and PTB controls (n=21); and (3) PRJEB30656 19 in which term FGR placental samples (n=5) were compared with term controls (n=5). For the GSE148241 data set, the VT blue, VT black, and DB brown modules were among those most strongly conserved, whereas the DB tomato and DB magenta modules were weakly preserved ( Figure S7C and S7D). The module preservation patterns for the GSE114691 samples were like those of GSE148241, except for the weak conservation of the VT Fountain module ( Figure S7C and S7D). The PRJEB30656 placentas showed evidence for retention of 4 VT modules (turquoise, black, straw, and blue) and 3 DB modules (mustard, magenta, and brown; Figure S7C and S7D).
Next, to gauge the utility of network modules in classifying clinical disease categories across studies (as may be applicable when detailed information on continuous clinical variables is inaccessible), we adapted and extended the metagene classification approach of Lauss et al. 9 A MGR was calculated for each module, representing the quotient of the averaged expression of the transcripts most positively and negatively correlated with specific continuous clinical variables ( Figures S8  and S9). These latter variables were selected based on their overall statistical association with each module in our specimens (Figures 2C and 3C).
Cluster analysis applied to the VT module MGRs revealed 3 major groups among the 4 data sets: a set with elevated VT black and VT blue MGR expression comprising a majority of the EOPE samples; a cluster containing a large proportion of control specimens; and an intermediate group ( Figure 4A). By ROC analysis, the VT black and VT blue MGRs showed the highest predictive performance for EOPE but had more limited associations with FGR ( Figure 4B and 4C). Although the VT orange module was significantly correlated with birthweight (r=0.60, Figure 2C) and had an MGR that was strongly associated with FGR in our samples (area under the ROC curve [AUC]=0.91, Figure 4C), its association with FGR in the other data sets was weak (AUC=0.61, Figure 4C). The VT turquoise MGR showed only limited predictive ability for EOPE ( Figure 4B) and lacked consistently significant associations with FGR ( Figure 4C).
Among DB module MGRs, unsupervised clustering demonstrated 4 major categories: a group comprising mostly control specimens having elevated DB tomato, DB mustard, and DB magenta MGRs with diminished DB brown MGR levels; a small cluster with reduced expression of all DB module MGRs; and 2 clusters containing a majority of the EOPE and FGR samples stratified by DB brown module MGR expression ( Figure S10A). By ROC analysis, the DB brown MGR was the most predictive for EOPE overall and exhibited some ability to detect FGR in the external data sets ( Figure S10B). Across studies, the MGRs of the DB magenta, DB mustard, and DB tomato modules exhibited significant associations with FGR in our data set but lacked consistent predictive performance more generally ( Figure S10C).

Targeted Pathway Analysis Unveils a Reproducible Transcriptional Signature in a Subset of EOPE Placentas Suggestive of a Glycolytic Metabolic Shift
In reviewing our results, we noted that cluster analysis of the VT black MGR transcripts subdivided the VT placental samples into distinct groups, including a set of EOPE±FGR samples with an expression pattern distinct from the remaining specimens ( Figures S8C  and S11A). Since exploratory functional enrichment analysis for this module showed an overrepresentation of gene sets involved in metabolism ( Figure 2F), we performed a targeted analysis of the genes in the glycolysis/gluconeogenesis (Kyoto Encyclopedia of Genes and Genomes [KEGG] hsa00010) pathway. These transcripts reproduced the clustering pattern observed for the VT black MGR transcripts, with a subgroup sharing an increased abundance of genes encoding key glycolytic enzymes (callout in Figure  S11B). Consistent with the observed enrichment for hypoxia-responsive pathways in the VT black module, this expression pattern suggested a bioenergetic shift away from aerobic oxidative phosphorylation towards increased lactate production. 22 To further simplify our analysis, we identified 2 RNA-seq expression ratios of VT black module genes: LDHA/LDHB; and PDK1/GOT1. The LDHA gene product preferentially converts pyruvate to lactate, whereas the LDHB gene product predominantly catalyzes the reverse reaction, 20 and the expression ratio of the 2 has been used as a surrogate index for oxidative phosphorylation (Warburg effect) usage in cancer cells. 21 The transcript for PDK1 is a direct target of HIF-1α and a central regulator of the pyruvate dehydrogenase complex, a gatekeeper in hypoxia-related metabolic reprogramming. 22 We prioritized GOT1 because it is repressed by HIF-1α experimentally 23 and displayed expression reciprocal to PDK1 in our data set (r=−0.76, P<0.0001). Clustering based on these ratios reiterated the distinct subgroup within EOPE±FGR samples seen with the KEGG pathway transcripts ( Figure  S11C). In the EOPE VT specimens, PC analysis of the VT black MGR and KEGG hsa0010 pathway transcripts revealed 2 distinct sample clusters ( Figure S11D) strongly correlated with the LDHA/LDHB and PDK1/GOT1 ratios ( Figure S11E). Division of the EOPE specimens into these 2 subgroups produced an expression pattern consistent with a glycolytic metabolic shift toward increased lactate production (cluster 1) and a pattern suggestive of normoxic glucose utilization (cluster 2; Figure S11F).The expression ratios in VT specimens by RNA-seq were significantly correlated with the corresponding ratios detected using qPCR (r=0.70, P<0.001 for LDHA/LDHB; r=0.56, P<0.002 for PDK1/GOT1; Figure S12). A comparison of the RNA-seq expression ratios within and between the VT and DB anatomical placental regions for the paired samples (n=27) is shown in Figure S13 where, in addition to the main VT cluster (cluster 1), a second cluster of 4 samples exhibiting selectively elevated DB ratios was also noted ( Figure S13D).
We broadened our analysis to the 3 previously published data sets to determine whether these findings were reproducible. Clustering of the batch-corrected, normalized expression data for the VT black MGR and KEGG hsa0010 transcripts consistently partitioned a common subset of samples (tan dendrogram lines, Figure 5A and 5B). This cluster could be reproduced using the RNAseq expression ratios LDHA/LDHB and PDK1/GOT1, and included the cluster 1 VT samples from our study ( Figure 5C). Across EOPE±FGR specimens from the available studies, combined analysis of these strongly correlated indices again revealed metabolic gene expression patterns consistent with a spectrum ranging from a standard bioenergetic profile to a shift similar to that seen in anaerobic glycolysis ( Figure 5E and 5F).

DISCUSSION
In this placental transcriptional profiling study, WGNCA revealed 14 correlated gene modules in our data set. All but 2 exhibited significant correlations with clinical variables (maternal blood pressure, birthweight, and GA at delivery) in our samples, providing a complementary alternative to differential abundance analysis for relating the placental transcriptome to EOPE and FGR. In considering the 160 samples from all 4 data sets, integrated analysis revealed 2 modules with strongly conserved EOPE and FGR associations: VT blue (1001 transcripts, highly correlated to BP and enriched for membrane receptor signaling pathways) and DB brown (124 transcripts, highly correlated to BP, comprised of transcripts involved in the regulation of the renin-angiotensin-aldosterone system, hypoxia-inducible factor, and steroidogenesis). These modules were 9.6% similar (by Jaccard index) and shared 99 genes in common, including several transcripts canonically associated with preeclampsia, such as LEP and FLT1. 24 Of the 5 modules most strongly correlated with birthweight, 4 (VT orange, DB tomato, DB magenta, and DB mustard) were significantly associated with FGR in our specimens but lacked strong associations when the 4 data sets were considered in aggregate. Our identification of the VT black module was of particular interest, given its independence from the larger VT blue cluster, its strong association with EOPE across studies, and its functional enrichment for transcripts associated with HIF signaling and core cellular metabolism; we, therefore, focused further consideration around this coexpression subnetwork.
Throughout pregnancy, the placenta must continually support the developing fetus while simultaneously executing its own energy-intensive biosynthetic activities; as such, it has a particularly high metabolic activity and requirement for oxygen. 25 Glucose is the main carbohydrate used by the uteroplacental system, and glycolysis serves as the central hub for placental metabolism. 25 Under stress conditions, particularly when oxygen supply is reduced, the placenta must compensatively reprogram its metabolism to keep pace with bioenergetic demands. In general, such circumstances require a switch to less energetically efficient pathways (eg, anaerobic glycolysis) than would be used otherwise. We thus anticipated that the placental insufficiency associated with EOPE and FGR-of which oxidative stress and uteroplacental hypoxia are cardinal features 6 -would involve the coordinated expression changes in the VT black module transcripts. However, we were surprised by the degree of variability in this expression signature among samples, both in our data set and in the EOPE specimens from public sources.
Closer examination of the glycolysis/gluconeogenesis pathway genes revealed a distinct molecular subgroup with a transcriptional signature consistent with a glycolytic shift; this cluster could be recapitulated using just the expression ratios of LDHA/LDHB and PDK1/GOT1 ( Figure 5). Strikingly, roughly half of all EOPE specimens (with or without FGR) and most of the isolated FGR specimens fell outside of this molecular cluster. To contextualize our present findings in light of the molecular classifications put forth by Leavey et al 11 and Gibbs et al, 12 we evaluated our data in comparison to the gene panels used by these investigators ( Figure S14). We approximated their canonical and immunologic subclasses by hierarchical clustering based on the relative expression of TAP1 to FLT1 and LIMCH1. Although an imperfect surrogate to their qPCR-based classification strategy, our analysis suggests that the metabolic shift signature occurred in both molecular subtypes described by Leavey et al 11 and Gibbs et al 12 and could represent an independent source of information within this framework.
It is now recognized that the plasticity associated with cancer cell metabolism is more nuanced than originally suggested by the classical model of Warburg-type glucose metabolism. 22 This is likewise the case with the intricacies of intermediary metabolism in the human placenta, 25 and the current results should be interpreted with these particularities in mind. It has recently been proposed that placental metabolic programming (eg, increased lactate production from glucose) may actually reflect a dynamic and adaptive physiological state, and that failure of such adaptations might contribute to the pathobiology of preeclampsia. 25 Although we cannot yet characterize our observations in terms of adaptation, decompensation, or malplacentation, it is possible that the altered glycolytic transcriptional patterns we identified may be transient, either spatially (through regional variation within the placenta) or temporally (eg, in response to acute stressors). Further work using animal models 26 and in vivo functional imaging of the human placenta 27 might help clarify the extent and nature of dynamic metabolic changes in complicated pregnancies.
Among the FGR cases, several coexpression modules were significantly associated with birthweight and FGR in our samples. However, these modules showed inconsistent evidence for preservation and poor associations with FGR in the other cohorts in which this condition was explicitly studied (GSE114691 and PRJEB30656). Indeed, we could identify no characteristic transcriptional signatures uniting these 3 studies. Research into FGR pathobiology is challenging given the diverse etiologies that might contribute to compromised fetal growth. 28 The clinical criteria used to diagnose FGR are continually evolving, and the distinction between constitutionally small for gestation age fetuses and those with FGR is not always clear. Therefore, it is unsurprising that prior placental transcriptional profiling studies of FGR have reflected this heterogeneity, particularly given that the criterion of small for gestation age alone has been frequently applied as a surrogate for FGR. 24 Although previous studies have reported overlap in the molecular processes between normotensive FGR and FGR associated with preeclampsia, 12,19,29 a recent large and comprehensive study reported multiple instances of divergence between FGR and EOPE expression patterns for several long RNA classes. 30 Since the isolated FGR cases were delivered preterm for fetal indications, it is impossible to predict whether the trajectory of these pregnancies might have included preeclampsia manifestations. Further, we cannot exclude the possibility that a small number of the FGR cases harbored associated (and potentially causal) chromosomal aberrations that were undetectable by traditional karyotyping. Further research to increase insight into the placental dysfunction associated with growth-compromised fetuses may require emphasis on pathway-level dysregulation rather than focusing on individual gene expression patterns.

Strengths and Limitations
Strengths of the present study include the availability of rich clinical phenotyping data for trait correlation analysis, integration with data from other studies to assess reproducibility, the separate consideration of VT and DB samples, and the use of GA-matched preterm controls in additional to term controls. The application of WGCNA as means for gene and preeclampsia sample classification has garnered increased research focus, [31][32][33] and it will be of interest to see how this approach may be further used for molecular subphenotyping of this syndrome in the near future.
Our study has several limitations. First, we recognize that isolated tissue biopsies may not reflect the global functioning of the placenta. We attempted to mitigate this potential shortfall by extending our observations to studies in which samples were pooled from multiple placental regions, finding that our principal observations were upheld in these other cohorts. Second, in our design, we focused on EOPE and FGR using carefully selected specimens that we believed would represent cases with the greatest potential to benefit from further targeted research; however, this restriction did not allow us to explore other clinical manifestations of these conditions. Finally, while powerful, we recognize that the analytical methods applied also entail shortcomings. 7 The foremost goal of this study was to explore the utility of correlation network analysis for prioritizing informative transcript modules in the settings of EOPE and FGR delivered before term. We recognize that identifying correlation networks more robustly would require a broader representation of the full clinical spectrum.

Perspectives
Our results contribute to a growing body of literature demonstrating that preeclampsia and FGR are more heterogeneous conditions than can be appreciated by applying clinical diagnostic criteria alone. We provide novel evidence for a molecular subphenotype consistent with a glycolytic metabolic shift that occurs frequently, but not universally, in placental specimens, suggesting a spectrum of placental responses to these clinical conditions. Translating this work into actionable methods to inform targeted therapeutics and guide personalized care will require further integration with key clinical variables including minimally invasive screening of placenta-derived materials, such as microparticles or cell-free nucleic acids. 34 Molecular subphenotyping through placental transcriptomics has promise to clarify and guide this work.