Developmental programming in human umbilical cord vein endothelial cells following fetal growth restriction

Background Fetal growth restriction (FGR) is associated with an increased susceptibility for various noncommunicable diseases in adulthood, including cardiovascular and renal disease. During FGR, reduced uteroplacental blood flow, oxygen and nutrient supply to the fetus are hypothesized to detrimentally influence cardiovascular and renal programming. This study examined whether developmental programming profiles, especially related to the cardiovascular and renal system, differ in human umbilical vein endothelial cells (HUVECs) collected from pregnancies complicated by placental insufficiency-induced FGR compared to normal growth pregnancies. Our approach, involving transcriptomic profiling by RNA-sequencing and gene set enrichment analysis focused on cardiovascular and renal gene sets and targeted DNA methylation assays, contributes to the identification of targets underlying long-term cardiovascular and renal diseases. Results Gene set enrichment analysis showed several downregulated gene sets, most of them involved in immune or inflammatory pathways or cell cycle pathways. seven of the 22 significantly upregulated gene sets related to kidney development and four gene sets involved with cardiovascular health and function were downregulated in FGR (n = 11) versus control (n = 8). Transcriptomic profiling by RNA-sequencing revealed downregulated expression of LGALS1, FPR3 and NRM and upregulation of lincRNA RP5-855F14.1 in FGR compared to controls. DNA methylation was similar for LGALS1 between study groups, but relative hypomethylation of FPR3 and hypermethylation of NRM were present in FGR, especially in male offspring. Absolute differences in methylation were, however, small. Conclusion This study showed upregulation of gene sets related to renal development in HUVECs collected from pregnancies complicated by FGR compared to control donors. The differentially expressed gene sets related to cardiovascular function and health might be in line with the downregulated expression of NRM and upregulated expression of lincRNA RP5-855F14.1 in FGR samples; NRM is involved in cardiac remodeling, and lincRNAs are correlated with cardiovascular diseases. Future studies should elucidate whether the downregulated LGALS1 and FPR3 expressions in FGR are angiogenesis-modulating regulators leading to placental insufficiency-induced FGR or whether the expression of these genes can be used as a biomarker for increased cardiovascular risk. Altered DNA methylation might partly underlie FPR3 and NRM differential gene expression differences in a sex-dependent manner.


Introduction
Fetal growth restriction (FGR) describes the condition in which the fetus fails to reach its genetically determined growth potential. FGR most commonly results from placental insufficiency, in which a reduced uteroplacental blood flow, oxygen, and nutrients toward the fetus lead to aberrant fetal growth. Compensatory physiological mechanisms are set in motion, such as fetal hemodynamic redistribution over organs and epigenetic alterations; while these fetal adaptations might improve fetal survival, they are considered to be unfavorable in the long run.
FGR is linked to an increased susceptibility for various noncommunicable diseases in adulthood, including cardiovascular and renal disease [1][2][3][4]. The Developmental Origins of Health and Disease (DOHaD) hypothesis states that this predisposition originates in the womb, when the adverse in utero environment influences epigenetic developmental programming [5][6][7]. Preclinical research strongly supports sex-specific programming of cardiovascular and renal disease in FGR offspring [8]. However, evidence for this concept has been less evident in humans [9,10]. Epigenetic differences have been observed in placental tissue and cord blood collected from pregnancies in which babies were born with a low birth weight (as a surrogate marker for FGR). Most studies investigated DNA methylation and identified epigenetic DNA methylation markers related to impaired fetal growth [11,12].
While placental tissue can be used to examine gene expression or epigenetic changes in pregnancies complicated by FGR, this tissue consists of a combination of maternal and fetal cells. Therefore, human umbilical vein endothelial cells (HUVECs) can be used to examine the fetal profile of disrupted growth without contamination by maternal cells. Therefore, HUVECs are especially relevant cell type in context of the fetal origin of cardiorenal disease. The few studies performed in cultured HUVECs report different proteome profiles in cultured HUVECs from FGR compared to control donors and differential protein expression and DNA methylation in the eNOS pathway [13][14][15]. However, transcriptomic profiling by RNA-sequencing of primary HUVECs derived from FGR compared to control donors, without culturing bias, has to our knowledge not yet been reported. In addition, the expression of gene sets related to cardiovascular or renal development and function has not been analyzed in this condition.
This study aims to examine whether developmental programming profiles, especially related to the cardiovascular and renal systems, differ in HUVECs collected from pregnancies complicated by placental insufficiencyinduced FGR compared to normal growth pregnancies. We explored this by whole-genome RNA-sequencing to map differential expression per gene and gene set enrichment analysis (GSEA) focussed on cardiovascular and renal development, function and health. Additionally, we performed targeted DNA methylation assays in candidate genes to gain insight in whether DNA methylation plays a regulatory role in the different expression. This approach contributes to the identification of early targets that can be aimed at to predict or prevent long-term diseases.

Study population
Pregnant women with placental insufficiency-induced FGR and pregnant woman with normal grown fetuses were included in this prospective observational study in the Wilhelmina Children's Hospital from July 2016 to December 2017. Inclusion criteria for placental insufficiency-induced FGR cases were described [16], but in short were diagnosed by prenatal ultrasound when (1) estimated fetal weight or abdominal circumference was below 10 th percentile for gestational age, in combination with (2) signs of placental insufficiency defined as abnormal blood flow patterns in umbilical artery, fetal middle cerebral artery, cerebral-placental ratio, or deflecting fetal growth rate in at least three consecutive measurements. The control group consisted of pregnancies with normal fetal growth defined as estimated fetal weight or abdominal circumference between p10-90. Percentiles of prenatal biometry were determined using the perinatology biometry calculator (http://www.perin atolo gy.com/ calcu lator s/biome try.htm). Exclusion criteria were congenital disorders, multiple pregnancies and stillbirth. The Medical Ethical Committee of the University Medical Center Utrecht approved the study on July 19, 2016, protocol number 16-302. Written informed consent was obtained from parents prior to delivery. preeclampsia, diabetes, autoimmune disorders. Use of maternal medication was registered. Percentiles for weight and head circumference at birth were determined with Intergrowth-21st [17]. Neonatal complications included infant respiratory distress syndrome, intraventricular hemorrhage, sepsis and necrotizing enterocolitis.

HUVECs isolation
Directly after placental delivery, the umbilical cord was stored in phosphate buffered saline (PBS) solution (pH 7.2;) at 4 °C. HUVECs isolation occurred preferably within 12 h, but always within 24 h after placental delivery as described [18]. Umbilical cords from n = 8 control and n = 12 FGR cases were collected. Cannulation of the umbilical vein at one end allowed access to wash with sterile PBS (pH 7.4; Gibco by Life Technologies, Grand Island, NY). Hereafter, the umbilical cord was clamped at both ends in order to incubate with accutase (0.02 µg/ml DNase; Innovative cell technologies Inc, San Diego, CA) for 5 min in sterile PBS at 37 °C to detach the endothelial cells from the umbilical vein. Detached HUVECs in accutase were flushed out of the umbilical vein with endothelial cell growth medium-2 (97% EGM-2; basal medium and SingleQuots supplement [1.9% FBS, 0.04% hydrocortisone, 0.4% hFGF-B2, 0.1% VEGF, 0.1% R3-IGF-1, 0.1% ascorbic acid, 0.1% hEGF, 0.1% GA-1000, 0.1% heparin], Lonza Bioscience, Walkersville, MD) and centrifuged in two separate tubes for 5 min 330 g at room temperature. One pellet was resuspended in 600 μl RA1 lysis buffer (Macherey-Nagel, Düren, Germany) and 6 μl 1 M DTT and stored at − 80 °C until RNA isolation. The second pellet was resuspended in 0.5 ml EMG-2 medium and 0.5 ml freezing medium with 20% DMSO and was frozen in a freezing container overnight and stored in liquid nitrogen the next day until DNA isolation.

RNA isolation and RNA-sequencing
RNA was isolated using NucleoSpin RNA ® (Macherey-Nagel), with RNA elution in 40 μl nuclease-free water. The concentration of RNA was quantified using Qubit RNA HS assay and Qubit fluorometer (Thermo Fisher). RNA-sequencing was performed as described [19]. In short, libraries were generated using NEXTFlexTM Rapid RNA-seq Kit (Bio Scientific) and sequenced by the Nextseq500 platform (Illumina) to produce 75 bp singleend reads through the Utrecht DNA sequencing facility. Reads were aligned to the human reference genome GRCh37 using STAR.

Gene set analysis
Gene set enrichment testing was performed on the hallmark (H), canonical pathway (C2-CP) and select GO term (C5) gene set collections from the Molecular Signatures Database (version 7.1) [20,21]. Only gene sets with relation to renal or cardiovascular development, function and health were selected from the GO term gene sets (Additional file 1: Table S1). Gene sets with less than five genes in the set of selected genes (based on expression, see below) were excluded from the analysis, eventually resulting in 2167 included gene sets.

DNA isolation and methylation
Genomic DNA from HUVEC was isolated with the allprep DNA/RNA mini kit (Qiagen, Venlo, the Netherlands), following the manufacturer's protocol. DNA quantity was measured with a Nanodrop 2000c (Thermo Scientific, Pittsburgh, PA). DNA was stored at − 80 °C until further analysis.
Targeted DNA methylation assaying was performed blindly in the significant differential expressed genes. Bisulfite conversion of 200 ng DNA was performed with the EZ DNA Methylation-Gold kit (Zymo Research, Leiden, the Netherlands) according to the manufacturers' protocol. Pyrosequencing primers were designed for the top three differentially expressed genes targeting the promoter regions (Additional file 2: Table S2)  The polymerase chain reaction product was analyzed for the extent of methylation per selected CpG positions on a PyroMark Q24 (Qiagen). Data were analyzed using the PyroMark Q24 Analysis Software 2.0 (Qiagen).

Statistical analysis Clinical data
Statistical analysis was performed in IBM SPSS Statistics 25 for Windows, version 25 (IBM Corp, Armonk, NY). Parametric data are presented as mean ± SD and tested with independent t test, nonparametric data are presented as median (minimum-maximum) and tested with Mann-Whitney, and nominal data are presented as n (%) and tested with Fisher exact. A two-sided p value of ≤ 0.05 was considered significant.

Differential expression of genes
Read counts per gene, per sample, were analyzed for global expression differences using R (version 3.5.3). Genes were selected with an expression of one count per million reads (CPM) in at least eight samples (n = 13,760 genes selected). Read counts were TMM-normalized using the calcNormFactors function from the edgeR package (version 3.24.3) [22]. TMM-normalized counts were used to assess global transcriptional profile differences of all samples by principal component analysis (PCA). Ten principal components (PC) were analyzed in the PCA analysis; values from each PC were checked for correlation to sample characteristics by the Mann-Whitney U test implemented in the SciPy package (version 0.19.0) in python (version 2.7.10). Outliers in RNA-seq data were identified and removed when (1) the number of reads was less than 1.000.000; (2) the number of nonzero genes was less than 10.000 or (3) a combination of number of nonzero genes was between 10.000 and 12.000 and a visible outlier on one of the PCA components. One outlier was identified; thus, n = 11 FGR and n = 8 control samples were selected for differential expression analysis.
Differential gene expression analysis was performed with the edgeR package (version 3.24.3) in R (version 3.5.3). Gene expression was modeled using the glmQLFit function in EdgeR [22], to a model that included patient group variables, as well as factors to capture Mode of Delivery (caesarean section vs. spontaneous delivery), Sex (male vs. female), and Gestational stage (preterm vs. term)-related gene expression variation. Differential gene expression was determined between study population groups (FGR vs. control). Differential expression statistics were obtained using the glmQLFTest functionality in edgeR; false discovery rates (FDRs) were determined using the Benjamini-Hochberg method to adjust for multiple testing and considered significant when below 0.1 (in combination with p value below 0.05) [23].

Gene set analysis
Gene set enrichment testing was performed with CAM-ERA, using the same linear model and contrasts as in the differential gene expression analysis (see above), and FDRs were also determined using the Benjamini-Hochberg method [23]. When a module showed ≥ 50% overlap with a higher ranking gene set, we selected the more significant gene set. Heatmaps for the gene sets related to the cardiovascular and renal development or function were created.

DNA methylation
The level of DNA methylation is given as a percentage, and since sex-specific differences have been reported in HUVECs and DNA methylation assays in other reproductive tissue, the data were analyzed with two-way ANOVA with Bonferroni multiple comparison using GraphPad Prism (version 8.4.3, San Diego, California, USA) [18,24,25].

Study characteristics
Study characteristics are presented in Table 1. There were no maternal cardiovascular diseases diagnosed besides preexisting hypertension or preeclampsia/HELLP. Severe FGR, defined as estimated fetal weight and/or abdominal circumference below the third percentile, was observed in ten out of 11 cases within the FGR group. One out of eight in the control group and four out of 11 in the FGR group were born prematurely. No neonatal death prior to discharged occurred. None of the neonates suffered from necrotizing enterocolitis or sepsis during neonatal intensive care unit admission, interventricular hemorrhage occurred in two control patients, and idiopathic respiratory distress syndrome was diagnosed in one control and one FGR neonate.

Table 1 Maternal and neonatal characteristics
Data expressed as mean ± SD or n (%), respectively, tested with independent t test or Fisher's exact test. # represents missing data, and therefore, the percentages are calculated based on the number of observations/ measurements within the control group with 7 being the lowest number of patients in a group (maximum 13% missing data). Pre-existing hypertension, preeclampsia and HELLP were defined according to the National Institute for Health and Care Excellence (NICE) guidelines [26]. GA: gestational age; HELLP: hemolysis, elevated liver enzymes and low platelet syndrome; and PPROM: preterm premature rupture of membranes

Differential expression of genes
Multidimensional scaling (MDS) plots showed clustering in the study population (FGR vs. CON) as potential modifiers, but not in delivery route, prematurity or sex (Additional file 3: Figure S1). PCA plots also showed clustering between FGR vs. CON (Additional file 4: Figure S2). All of the study characteristics were tested for association for all the first ten PCs, and study population was associated with PC1, PC3, PC4 and PC6; delivery route was associated with PC1 and PC7, and gestational stage with PC6 (Additional file 5: Table S3). Therefore, differences in expression due to sex, mode of delivery, and gestational stage were accounted for in the modeling of gene expression. Three protein-coding genes and one long intergenic noncoding (linc)RNA gene were significantly regulated (with a FDR < 0.1; Additional file 6: Table S4) in FGR compared with control samples: 1) lectin, galactosidebinding, soluble, 1 (LGALS1), 2) formyl peptide receptor 3 (FPR3), 3) nurim nuclear envelope membrane protein (NRM), 4) lincRNA RP5-855F14.1; all protein-coding genes were downregulated and lincRNA gene was upregulated ( Fig. 1).

Differential expression of gene sets
In total, 336 gene sets were significantly different between FGR and CON (with a FDR < 0.1; Additional file 7: Table S5). Selection of only the highest ranking gene set module (overlapping modules excluded) resulted in 193 downregulated and 22 upregulated gene sets in FGR versus control samples. The downregulated gene sets are mostly involved with immune, inflammatory or cell cycle pathways. As we were interested in risk of developing cardiovascular or renal disease, we noticed that the several gene sets related to renal development were significantly upregulated and a few related to cardiovascular health and function were downregulated in FGR samples vs. CON samples ( Table 2). Heatmaps were made to study the extent of up and downregulation for the distinct genes in these gene sets in each sample (Additional file 8: Figure S3). From this analysis, most genes were up-and downregulated in accordance with the differential gene set analysis results.

DNA methylation
The percentage methylation was measured at each individual CpG position for selected areas of the promoters of the three protein-encoding genes LGALS1, FPR3 and NRM per study population (Table 3).
LGALS1 showed similar percentage of methylation at each individual CpG position between groups, independent of sex (Fig. 2). DNA hypomethylation differed between in FGR males versus CON males at CpG2 of FPR3 only (Additional file 9: Figure S4). NRM was significantly hypermethylated at CpG1 in FGR compared to control, especially in male FGR offspring (Fig. 3).

Discussion
This study examined developmental cardiovascular and renal programming profiles in HUVECs collected from pregnancies complicated by placental insufficiencyinduced FGR compared to normal growth pregnancies to identify targets underlying long-term cardiovascular and renal diseases. We report downregulated expression of the protein-coding genes LGALS1, FPR3 and NRM and upregulation of the lincRNA RP5-855F14.1 in FGR samples compared to controls. Sex-dependent DNA methylation might partially underlie FPR3 and NRM gene expression, but we did not observe this for LGALS1. Additionally, of the significantly differentially expressed gene sets, the downregulated ones were mostly involved with immune, inflammatory, or cell cycle pathways. Interestingly, seven of the 22 significantly upregulated gene sets related to kidney development and four gene sets related to cardiovascular function and health differed between FGR and control.

Downregulated expression of LGALS1
LGALS1 is the protein coding gene for galectin-1 (Gal-1) [27]. During pregnancy, Gal-1 is important for immunomodulatory and vascular adaptions required for healthy placentation [28][29][30]. In HUVECs, Gal-1 mediates angiogenesis via vascular endothelial growth factor receptor (VEGFR)-2 but also the neuropilin receptor (NPR)-1 which enhances the binding between VEGF and LGALS1 expression was downregulated in placentas derived from early-onset preeclamptic patients, which is a placental insufficiency syndrome just as in FGR and part of our study population [30]. In addition, in a mouse model inhibition of Gal-1 mediated angiogenesis resulted in preeclamptic symptoms and fetal growth restriction [30]. However, higher or no difference in Gal-1 expression was observed in term placentas from pregnancies complicated by respectively preeclampsia or FGR [32,33]. Considering the different histopathology between early-onset and late-onset FGR (or preeclampsia), the downregulated LGALS1 expression in early-onset FGR speculatively contributes directly to placental insufficiency, while the upregulated expression observed in late-onset FGR might be secondary to relative placental-umbilical hypoxia or reduced umbilical flow [30,32,34]. Gal-1 has been suggested as an early marker of endothelial dysfunction, and dysregulated Gal-1 has been linked to poor blood pressure regulation and development of cardiovascular disease [27,28]. Therefore, the finding of downregulated LGALS1 expression in our FGR samples might be the key regulator leading to placental insufficiency-induced FGR, as well as an indication for the possible higher risk of developing long-term cardiovascular dysfunction in these offspring. Epigenetic processes or post-transcriptional modifications other than DNA methylation might be involved in reduced LGALS1 expression.

Downregulated expression of FPR3 and NRM and upregulation of lincRNA RP5-855F14.1
In contrast to LGALS1, FPR3, NRM, and RP5-855F14.1 are less studied genes, especially in the context of FGR and pregnancy. To our knowledge, this is the first study  reporting downregulated expression of FPR3 and NRM and upregulation of RP5-8551F14.1 in HUVECs or other pregnancy tissues collected from pregnancies complicated by FGR. One study in mice suggested that NRM might be involved in early cardiac development [35]. While lincRNAs in general have been described to be epigenetic process-associated with cardiovascular disease and development, the exact function of this specific lincRNA is unknown [24,36]. In humans, circulating FPR3 mRNA in combination with three other circulating mRNA had a high specificity and sensitivity to predict early-onset PE [37], but little is known on FPR3 expression or function. Recent studies did link downregulated expression of FPR2-another member of the FPR family with high analogies in sequencing identity and downstream responses [38]-in placental tissue to endothelial dysfunction and placental insufficiency via impaired immunomodulatory and angiogenic processes, leading to FGR [39,40]. FPR2 has also been described to play a protective and repairing role in ischemic heart disease and stroke [41,42]. Considering the similarities between FPR2 and FPR3 [38], the observed downregulation of FPR3 in our placental insufficiency-induced HUVECs samples could potentially contribute to placental insufficiency in a similar manner and dysregulated FPR3 might contribute to the increased susceptibility to cardiovascular disease. DNA methylation might partly contribute to the differential gene expression patterns of FPR3 and NRM in a sex-dependent manner. The sex dependency is especially interesting given that hypertension is more pronounced in fetal growth restricted male offspring [6,8,43]. However, while the methylation differences are significant for NRM, they are relatively small suggesting that the regulating role of DNA methylation in NRM expression is not strong. In addition, FPR3 is expressed, while the CpGs show hypermethylation; since the closest CpG position is almost 400 bp upstream of the TSS (although not by exclusion), its methylation possibly has little regulatory effects.

Cardiovascular and renal gene sets
We focused our gene set analysis on cardiovascular and renal development and function, since FGR has been associated with increased susceptibility to develop cardiovascular and renal disease. The four different cardiovascular gene sets, including lipid metabolism, might be in line with the cardiovascular risk profile described with dysregulated expression of LGALS1 and RP5-8551F14.1. While the gene sets related to cardiovascular development were similar in both groups, several gene sets related to kidney development were relatively upregulated in FGR. Nephrogenesis starts around 22 days of gestation and is complete at 34-36 weeks of gestation, making pregnancy the most vulnerable period to impact nephron endowment [44]. FGR has been linked to reduced nephron count and morphological differences in glomeruli, which could lead to glomerular hypertension and compensatory hyperfiltration in the remaining nephrons, causing subsequent nephron loss (Brenner hypothesis) [45]. A recent study using three different rat models for FGR also reported that molecular pathways differed in kidneys from FGR and control male offspring at birth and at postnatal day seven (end of nephrogenesis in rats) [46]. The pathways involved depended on the stage of development, and most upregulation was observed in the placental insufficiency-induced model (best matching our study population). The upregulation of the gene sets related to kidney development in our study, although not evident at individual gene level, suggests that developmental programming difference as a consequence of FGR is a subtle process. In an animal study, upregulation of renal genes in combination with wider nephrogenic zones suggested delayed nephrogenesis in FGR [47]; however, due to ethical reasons, we cannot confirm this histologically. Therefore, whether the upregulation in HUVECs represents accelerated or delayed kidney development and whether this relates to long-term renal function or disease remain to be elucidated.

Strengths and limitations
This is the first study using HUVECs to investigate developmental programming differences between placental insufficiency-induced FGR and controls by a full-transcriptome RNA-sequencing approach followed by differential gene and gene set expression analysis and follow-up DNA methylation assaying. The gene set analysis focused on cardiovascular and renal development, to specifically test the hypothesis that this correlate with increased risk for diseases of these systems can already be found during pregnancy. The transcriptomic profiling by RNA-sequencing enabled us to find FGR associated gene expression differences in an unsupervised manner, and to select genes for the analysis of DNA methylation. A major strength is that we used prenatal ultrasound measurements to clearly define the placental insufficiency-induced FGR phenotype in our study population; most studies use birth weight as surrogate marker, but this umbrella term also includes other underlying mechanisms such as congenital disorders or constitutionally small children, which have not been exposed to a hostile in utero environment. The strength of investigating HUVECs is that they allow us to investigate fetal expression differences without the interference of maternal cells. Additional strengths are that we did not culture our collected HUVECs and that we isolated the HUVECs within 24 h after delivery which limits external influences on epigenetic results. While our strictly defined study population created a clear placental insufficiency-induced FGR phenotype, it also limited the number of included samples. However, (post hoc) power calculations for RNA-seq data are not well developed and as such not custom to use for this type of analysis. Instead, we relied on multiple testing corrects, and variation estimates build into the edgeR pipeline to properly analyze our data. Although beyond the scope of our study, we acknowledge that we cannot correlate gene expression findings to long-term outcomes. Gestational age at birth differed between groups, but we accounted for this factor in modeling gene expression and plots of gene expression per sample versus gestational age at birth showed that this was not of influence. The downside of not culturing our HUVECs was that the yielded RNA concentrations were relatively low. However, careful consideration of the number of reads and nonzero genes in all samples allowed detection and exclusion of low-quality samples. Despite extensive washing, HUVEC samples might have been contaminated with a few other fetal blood cells. Other (epi)genetic mechanisms besides DNA methylation possibly involved in differentially regulated expression of genes were not tested. Confirmation of our programming hypothesis includes verification in other tissues.

Summary and future perspectives
In conclusion, this study showed downregulated expression of LGALS1, FPR3 and NRM and upregulation of RP5-855F14.1 in HUVECs collected from placental insufficiency-induced FGR compared to control. Additionally, several gene sets related to kidney development were upregulated and a few gene sets related to cardiovascular risk were downregulated in FGR. How these findings correlate with long-term cardiovascular and renal function requires further investigation and follow-up studies. The differentially expressed genes (or their encoded protein) might be used as a biomarker, which could contribute to personalized care by predicting the risk of developing cardiorenal disease and selective follow-up of only the patients at risk. Further studies are also required to elucidate how and whether the downregulation of LGALS1 and FPR3 are causal regulators resulting in placental insufficiency-induced early-onset FGR as they might hold promise as potential novel targets for preventive treatment. While the findings of this study are the first to support the concept of developmental epigenetic programming of cardiorenal disease following placental insufficiency-induced FGR in humans, they are also merely the first of many steps toward clinical applicability.