Abstract
In the human placenta, two trophoblast cell layers separate the maternal blood from the villous basement membrane and fetal capillary endothelial cells. The inner layer, which is complete early in pregnancy and later becomes discontinuous, comprises the proliferative mononuclear cytotrophoblasts, which fuse together and differentiate to form the outer layer of multinucleated syncytiotrophoblasts. Because the syncytiotrophoblasts are responsible for key maternal-fetal exchange functions, tight regulation of this differentiation process is critical for the proper development and the functional role of the placenta. The molecular mechanisms regulating the fusion and differentiation of trophoblasts during human pregnancy remain poorly understood. To decipher the interactions of non-coding RNAs (ncRNAs) in this process, we exposed cultured primary human trophoblasts to standard in vitro differentiation conditions or to conditions known to hinder this differentiation process, namely exposure to hypoxia (O2 < 1%) or to the addition of dimethyl sulfoxide (DMSO, 1.5%) to the culture medium. Using next generation sequencing technology, we analyzed the differential expression of trophoblastic lncRNAs, miRNAs, and mRNAs that are concordantly modulated by both hypoxia and DMSO. Additionally, we developed a model to construct a lncRNA-miRNA-mRNA co-expression network and inferred the functions of lncRNAs and miRNAs via indirect gene ontology analysis. This study improves our knowledge of the interactions between ncRNAs and mRNAs during trophoblast differentiation and identifies key biological processes that may be impaired in common gestational diseases, such as fetal growth restriction or preeclampsia.
Introduction
The fusion of mononucleated cytotrophoblasts into multinucleated syncytiotrophoblasts is a central process in human trophoblast differentiation. Early in pregnancy, this fusion process is a part of the pre-lacunar and lacunar stages of implantation on days 6–12 after fertilization in human pregnancy (). Once villi are formed, the fusion of mononucleated cytotrophoblasts into overlying multinucleated syncytiotrophoblasts at the villous surface is accompanied by a dramatic change in cell morphology, transcriptional output, and the production of growth factors and endocrine signals (). Located at the surface of human placental villi, the syncytiotrophoblasts are uniquely positioned to regulate key functions of the placenta in terms of maternal-fetal gas exchange, the uptake of nutrients into the feto-placental compartment, the release of waste to the maternal blood, the production of hormones and the immune and mechanical protection of the developing fetus (; ). The subjacent mononucleated cytotrophoblasts, which form a continuous layer early in pregnancy, later become a discontinuous layer of interspersed cytotrophoblasts that function as progenitors for replenishment of damaged or dead syncytium and homeostatic preservation of this critical layer (). The syncytiotrophoblast exhibits polarity, with a microvillous plasma membrane facing the maternal blood on the apical side and a basal plasma membrane located adjacent to the cytotrophoblasts and the basement membrane. Considering its functions, it is not surprising that the syncytiotrophoblast microvillous membrane, which interfaces directly with the maternal blood, harbors receptors for diverse plasma proteins, growth factors, immunoglobulins, and other soluble ligands, all linked to intracellular trophoblast signaling cascades (; ).
Syncytium formation prior to 10–13 weeks of human pregnancy takes place in a hypoxic environment (; ; ). Beyond that point, and once the maternal blood begins to perfuse the intervillous space, hypoxia may be harmful for proper placental function, leading to cytotrophoblast proliferation, attenuated fusion of cytotrophoblasts into syncytiotrophoblasts, reduced hormone and other biosynthetic functions, and overall syncytial damage and trophoblast death, leading to placental injury and diseases such as fetal growth restriction (; ; ; ; ; ; ; ; ; ). Seeking to characterize gene expression changes that define trophoblast differentiation, researchers have focused on the effect of hypoxia on the expression of protein-coding genes in term trophoblasts (; ; ; ).
Recent progress in untangling the complexity of the RNA world has shed light on diverse non-coding RNAs (ncRNAs) that play an essential role in shaping cellular differentiated functions. Among these RNAs, microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) represent the two best characterized classes and perhaps those with the most important regulatory potential. Over 2000 miRNAs are encoded in the human genome. Most of these miRNAs act in the cytosol, where they target mRNAs though imperfect base-pairing to block their translation and accelerate their decay (). However, despite their relative simplicity, the full impact of miRNAs on gene expression remains incompletely understood. lncRNAs are more diverse, with an estimated 30,000–100,000 nuclear and cytoplasmic species expressed from the human genome (; ; ; ). Further, the action of lncRNAs is complex, spanning interactions with DNA, RNA, and proteins and involving 3D structural flexibility that enables protein scaffolding and the assembly of multi-subunit complexes and nuclear condensates that shape transcriptional and posttranscriptional functions (; ).
Recent discoveries within the field of placental biology highlighted the putative role of lncRNAs and miRNAs in trophoblastic gene regulatory networks, their role in trophoblast differentiation and in response to hypoxic injury, and the impact of these processes on clinically relevant placental diseases (; ; ). A systematic inquiry into network interactions of lncRNAs, miRNAs, and mRNAs in differentiating primary human trophoblasts (PHT cells) is lacking. Here, we used an in vitro model of cultured PHT cells to investigate harmonized changes of lncRNAs, miRNAs, and mRNAs during PHT cell differentiation. In addition to exposure of PHT cells to hypoxia (), diverse chemicals and culture conditions have been employed in vitro to modulate the differentiation of PHT cells. Douglas et al. found that cytotrophoblasts exposed to 1.5% DMSO retain their mononuclear morphology, culminating in drastic inhibition of hCG production (). Other approaches to limit trophoblast differentiation include the use of colchicine, an inhibitor of microtubule polymerization (), cobalt chloride, a hypoxia-mimicking agent (; ), and the use of Ham’s/Waymouth medium (; ; ). The deployment of these culture conditions led to the discovery of a repertoire of genes that have been implicated in trophoblast differentiation and to which hypoxia-induced placental injury has been attributed (; ; ).
We sought to hinder trophoblast differentiation using hypoxia or the addition of DMSO to the culture medium, two approaches that have led to reproducible results in our laboratory (; ; ; ; ; ; ). We used next generation sequencing technology to identify differentially expressed lncRNAs, miRNAs, and mRNAs during these processes. Importantly, we interrogated the interactions among these RNAs and inferred the main biological processes represented by co-expression patterns.
Materials and Methods
Placentas and Dispersed Primary Human Trophoblasts
All placentas used in our studies were obtained from uncomplicated pregnancies and term deliveries at Magee-Womens Hospital in Pittsburgh, under a protocol that was approved by the Institutional Review Board at the University of Pittsburgh. PHT cells were isolated using the trypsin-DNase-dispase/Percollx method as described by , with modifications as we previously published (; ). PHT cells were cultured in DMEM (Sigma-Aldrich, St. Louis, MO) containing 10% bovine growth serum (HyClone, Logan, UT) and 1% P/S antibiotics (Sigma-Aldrich) at 37°C in a 5% CO2-air atmosphere, until culture conditions were modified as below.
The data used for this study were derived from two independent sets of experiments, each including several paradigms (Figure 1). In the first experimental set, PHT cells from five independently collected placentas were first cultured for 4–6 h in standard conditions (20% O2) to allow adhesion, using protocols established in our lab. Some of the cells were harvested at the end of the initial incubation period (time 0 in Figure 1) and were used as control. The remaining plates from the same culture were maintained for an additional 48 h in either standard culture conditions or in hypoxia (O2 < 1%), using a dedicated hypoxia chamber, as we previously described (). In the second experimental set, PHT cells from six independently collected placentas were first cultured for 4–6 h in standard conditions as above. The culture continued for an additional 48 h, with some of the cells exposed to DMSO 1.5% (Sigma-Aldrich), designed to mitigate cell differentiation as previously shown () and reproducibly validated by us (; ).
FIGURE 1
RNA Extraction, Library Preparation, and Sequencing
Total RNA was isolated from from placental specimens by using TRI Reagent (Sigma) according to the manufacturer’s instructions and purified using EconoSpin spin columns (Epoch Life Science, Missouri City, TX). The quantity and quality of total RNA was determined with a NanoDrop 1000 spectrometer (Thermo Fisher, Waltham, MA) and an Agilent bioanalyzer (Agilent Technologies, Santa Clara, CA). From each extracted RNA sample, we used 10 μg of total RNA to generate two types of libraries—one for long RNAs (≥200 nt), including mRNAs and lncRNAs, and one for small RNAs. The libraries were prepared and sequenced by Ocean Ridge Biosciences (Palm Beach Gardens, FL) and by McGill University’s Génome Québec Innovation Centre (Montréal, Canada). For the second set of experiments, the libraries were prepared and sequenced at the Health Sciences Sequencing Core at Children’s Hospital of Pittsburgh. The miRNA samples were sequenced using the QIAseq miRNA sequencing protocol, which links a unique molecular identifier (UMI) to each miRNA to reduce sequencing bias. Data from all experiments were deposited to the Sequence Read Archive (SRA) at the National Center for Biotechnology Information with BioProject IDs: PRJNA674312, PRJNA674329, PRJNA674366, PRJNA704383, PRJNA704399, PRJNA704393.
Reverse Transcriptase and Quantitative PCR (RT-qPCR)
RNA was extracted from cells with TRI Reagent (Sigma). cDNA was synthesized from 1 μg of total RNA by using the High-Capacity cDNA Reverse Transcription kit (Applied Biosystems, Foster City, CA) according to the manufacturer’s protocol. Template cDNA was PCR-amplified with the gene-specific primer sets (Supplementary Figure 1B). RT-qPCR was performed using SYBR Select (Applied Biosystems) in a ViiA 7 system (Applied Biosystems). Analysis of qPCR data was performed using the delta-delta Ct method (), normalized to glyceraldehyde-3-phosphate dehydrogenase expression.
RNAseq Data Processing
The long RNA libraries were aligned to human reference genome GRCh38 using STAR (2.5.2b), an RNAseq alignment tool (), and annotated with GENCODE v25 (). The number of reads per gene was calculated for each RNAseq library, using STAR. We used the definition of lncRNA, which includes antisense RNA, sense intronic RNA, processed transcripts, and sense-overlapping RNAs (). The 15 small RNA libraries from the first set of experiments were analyzed using our validated miRNA sequencing data analysis pipeline (). Briefly, after pre-processing the library reads, including the removal of adaptor sequences and dimerized primer sequences, we used Bowtie to align all remaining reads with of least 15 nucleotides to the human reference genome (GRCh38) (). The BEDTools program was used to intersect the alignments to a mature miRNA database maintained by miRBase (v21). The intersected alignments were summarized to obtain counts for all miRNAs. The 18 small RNA libraries from the second set of experiments were processed using the online Qiagen Primary QIAseq miRNA quantification tool1 to handle the reads with UMIs. We found an extremely high correlation between miRNA counts with distinct UMIs and total miRNA reads. To ensure consistency across the two experimental sets, the total miRNA reads were used as the miRNA expression data.
The counts of long and short RNAs in the sequencing libraries were assumed to follow negative binomial distributions. The negative binomial test, implemented in the Bioconductor R package DESeq2 (), was used to identify differentially expressed lncRNAs, miRNAs, and mRNAs. Fisher exact test was used to perform gene ontology analysis and identify the biological processes that were over-represented in selected mRNAs. All p-values from multiple simultaneous multiples tests were adjusted using Benjamin and Hochberg’s method to control the false discovery rate ().
Model-Based Co-expression Analysis
We proposed a model-based method to identify mRNAs that were co-expressed with lncRNAs and/or miRNAs. Specifically, we assume that the two genes are co-expressed if their expressions Y and X are related through a (generalized) linear model:
where g is the link function, Z is the vector of confounding variables, B is a row vector representing the coefficients of Z, and a is a non-zero coefficient representing the co-expression relation between X and Y. We then examined whether X and Y were co-expressed by testing if a = 0. Note that when X, Y, and Z have a joint multivariate normal distribution, this is equivalent to testing whether X and Y have a zero partial correlation, given Z. As we used a negative binomial model to analyze the gene expression, Y is the count for one gene and X the log-transformed and normalized expression of the other gene. The vector B accounts for the experimental exposure and for placentas used for each batch of PHT cells.
Analysis of Gene Ontology
Because the functions of lncRNAs and miRNAs are not available from the gene ontology database, we proposed an indirect approach through the gene co-expression analysis described above. Specifically, given a set of lncRNAs or miRNAs and using the model-based co-expression analysis, we identified all mRNAs that are co-expressed with each of the lncRNAs or miRNAs. The top-ranked mRNAs were then used for analysis of gene ontology to identify the enriched biological processes, attributed to altered expression of the set of lncRNAs or miRNAs.
The statistical analyses, detailed above, were performed using R () and Bioconductor (). For analysis of lncRNA by RT-qPCR, the fold-change data were analyzed using Kruskal Wallis non-parametric test, with post hoc Tukey test for all pairwise comparisons. RT-qPCR data analysis was performed using Prism software (GraphPad Software, San Diego, CA).
Results
Using an average read of 0.5 per library as a cutoff in the first set of experiments, we identified 33,719 long RNAs, including 8,035 lncRNA and 17,298 mRNA species, from the 15 long RNA libraries and 923 miRNAs from the 15 small RNA libraries. In the second set of experiments, we identified 30,421 long RNAs from the 15 long RNA libraries, including 6,988 lncRNA and 16,723 mRNA species, and 2,416 miRNAs from the 15 small RNA libraries. Overall, mRNAs had a higher read count per library than lncRNAs (Figure 2). For example, in the first experiment the median of average lncRNA reads per library was 6.2, compared to 457.5 reads for the mRNAs. In the second experiment, the corresponding values were 4.2 and 267.7, respectively. We also used RT-qPCR to validate the expression changes of 10 lncRNAs (Supplementary Figure 1), as we previously did for mRNA and miRNA transcript data (; ).
FIGURE 2
To visualize the effect of exposure and potential confounding factors on the PHT transcriptome, we normalized and log transformed the long RNA and miRNA library data, using the regularized and variance-stabilizing transformation method in the Bioconductor package DESeq2. We then performed classical multidimensional scaling, respectively, for lncRNAs, miRNA, and mRNAs (Figure 3). These plots showed that the experimental conditions markedly influenced the PHT transcriptome. The exposure effect dominated for mRNAs in the two experimental sets and for lncRNAs in the second set and was clearly visible for lncRNAs in the first set and for miRNAs in both experimental sets. The plots also showed a batch effect on the expression of lncRNAs and miRNAs in the first experimental set (hypoxia), as samples 1 and 2 were processed by a different lab than samples 3–5. Our data also suggest an effect of the placenta on miRNAs in both experimental sets.
FIGURE 3
We tested the lncRNAs, miRNAs, and mRNAs that were differentially expressed between the conditions in the two experimental sets. For this, we used negative binomial regression models with two factors: the exposure factor (hypoxia, DMSO) and the placenta factor, which refers to inter-individual variation in placental transcriptome and possible differences in processes related to placental collection and cell preparation. After controlling the false discovery rate at 0.05, the numbers of differentially expressed lncRNAs, miRNAs, and mRNAs are shown in Table 1. For convenience, herein genes differentially expressed under the hypoxic as opposed to the standard condition are termed “hypoxia-DE” genes, and genes differentially expressed with DMSO as distinguished from the standard condition are termed as “DMSO-DE” genes. We next examined the number of lncRNAs, miRNAs, and mRNA that exhibited concordant (Figure 4) or discordant (Supplementary Figure 2) expression change across the two experimental sets. We noticed a significant concordance (up or down) in RNA expression change between the two sets of experimental exposures, namely “hypoxia-DE” and “DMSO-DE” genes (Figure 5). This concordance was noted primarily for lncRNA and mRNA species (Figures 5A,C), with less concordance for up- or downregulated miRNAs (Figure 5B).
TABLE 1
| Std vs. Hypoxia | Std vs. DMSO | |||
| Up | Down | Up | Down | |
| mRNA | 4,803 | 4,236 | 4,354 | 4,607 |
| lncRNA | 571 | 1,658 | 1,012 | 422 |
| Other RNA | 311 | 1,234 | 520 | 264 |
| miRNA | 69 | 58 | 94 | 156 |
The numbers of differentially expressed RNA transcripts in the two experimental sets.
FIGURE 4
FIGURE 5
We also assessed the similarity of the two experimental sets by examining expression differences in mRNAs, lncRNAs, and miRNAs, between the 0 and 48 h time points, under standard conditions. We found an extremely high concordance (up or down) in the expression change of all three types of RNAs between the two experimental sets (Supplementary Figure 3). This confirms that the two experimental sets are comparable for similar conditions, thus supporting our experimental approach.
We provided additional support to our findings by comparing our gene expression changes with published human placental single cell RNAseq (scRNAseq) data, predicting that our data would be similar to gene expression changes between syncytiotrophoblasts and cytotrophoblasts. We performed this analysis using the PlacentaCellEnrich tool (), available for protein coding genes, and examined changes among the 19 mRNAs with the highest concordant log2-fold increased expression and the 6 mRNAs with the highest concordant log2-fold decreased expression between standard conditions and hypoxia/DMSO. Among the selected 25 mRNAs (Supplementary Table 1) the vast majority of our results were consistent with the scRNAseq data, as predicted, with expression changes between standard conditions vs. hypoxia/DMSO that are similar to changes between syncytiotrophoblasts and cytotrophoblasts. Notably, discrepancies were more common for one of the databases, likely reflecting much lower mRNA expression levels and a lower magnitude of the log2-fold change for mRNAs that are downregulated during cell differentiation.
To define pathways that might be implicated in reduced trophoblast differentiation, we performed direct and indirect gene ontology analyses. We focused on genes with at least a moderate expression, defined as having a minimum of 25 reads per library. For mRNA, we focused on 3,338 hypoxia-DE genes, 2,232 DMSO-DE genes, and 988 DE genes that were concordantly altered in the two sets of experiments. Table 2 shows the shared biological processes that were significantly enriched based on mRNA changes in “hypoxia-DE” and in “DMSO-DE” transcripts, irrespective of lncRNA or miRNA co-expression. Supplementary Table 2 shows differentially expressed mRNAs that were concordantly expressed in the two exposure sets and that are co-expressed with at least one differentially expressed lncRNA and miRNA, where both lncRNA and miRNA exhibit concordant expression change (up- or downregulation) between the two exposure sets.
TABLE 2
| Process |
| Immune response |
| Female pregnancy |
| Signal transduction |
| Inflammatory response |
| Cell surface receptor signaling pathway |
| Immune system process |
| Chemokine-mediated signaling pathway |
| Cell-cell signaling |
| Response to drug |
| Retinoid metabolic process |
| Cytokine-mediated signaling pathway |
| Positive regulation of cytosolic calcium ion concentration |
| Positive regulation of T cell proliferation |
| Defense response |
| Regulation of calcium ion-dependent exocytosis |
| Calcium ion-regulated exocytosis of neurotransmitter |
| Multicellular organism development |
| Cell adhesion |
| Positive regulation of phagocytosis, engulfment |
| Chemotaxis |
| Single organismal cell-cell adhesion |
| Innate immune response |
| Keratinization |
| Cellular response to lipopolysaccharide |
| G-protein coupled receptor signaling pathway |
| Positive regulation of ERK1 and ERK2 cascade |
| Positive regulation of cell proliferation |
Biological processes enriched in mRNAs that exhibited a concordant expression change during trophoblast differentiationa and ranked by adjusted p-value.
aTrophoblast differentiation was defined as expression change in control vs. hypoxia and in control vs. DMSO exposure.
To further pursue the biological pathways reflecting mRNA expression changes that were associated with differentially expressed lncRNAs and/or differentially expressed miRNAs, we performed indirect gene ontology analysis as described in section “Materials and Methods.” We first selected the 438 hypoxia-DE lncRNAs and 85 hypoxia-DE miRNAs, 167 DMSO-DE lncRNAs and 115 DMSO-DE miRNAs, as well as 66 concordant DE lncRNAs and 9 concordant DE miRNAs. We then performed model-based co-expression analysis for these associations on the basis of our negative binomial regression model. Table 3 shows biological processes attributed to mRNAs that were expressed in a concordant manner in the two exposures, in association with changes with either lncRNA or with miRNA in the same conditions. Lastly, we identified the biological processes attributed to mRNAs that were differentially expressed in a concordant manner between the two exposures in association with changes with both lncRNA and with miRNA in the same conditions (Table 4). These associations are also shown in the Circos plot (Figure 6).
TABLE 3
| Association | Biological processes |
| mRNA associated with lncRNAa | Positive regulation of transcription from RNA pol II promoter |
| Response to drug | |
| Inflammatory response | |
| Response to mechanical stimulus | |
| Positive regulation of gene expression | |
| Cellular oxidant detoxification | |
| mRNA associated with miRNAa | Female pregnancy |
| Cell adhesion | |
| Inflammatory response | |
| Extracellular matrix disassembly | |
| Cell-cell signaling | |
| Positive regulation of angiogenesis | |
| Response to drug |
Biological processes enriched in mRNAs that were associated with lncRNA or miRNA during trophoblast differentiation, and ranked by adjusted p-value.
aAt least one DE lncRNA or one DE miRNA, respectively, were associated with mRNA in the two exposures.
TABLE 4
| Process |
| Female pregnancy |
| Cholesterol metabolic process |
| Lipoprotein metabolic process |
| Positive regulation of cytosolic calcium ion concentration |
| Placenta development |
| Cell adhesion |
| Inflammatory response |
| Positive regulation of phagocytosis, engulfment |
| Retinoid metabolic process |
| Cytokine-mediated signaling pathway |
| Cellular response to hormone stimulus |
| Signal transduction |
| Immune response |
| Chemotaxis |
| Response to drug |
Biological processes enriched in mRNAs that were associated with lncRNA and miRNA during trophoblast differentiationa, and ranked by adjusted p-value.
aAt least one DE lncRNA and one DE miRNA were associated with altered mRNA expression in the two experimental exposures.
FIGURE 6
Discussion
To decipher RNA interactions during differentiation of primary term human trophoblasts, we interrogated two sets of experiments, where the normal differentiation process was hindered either by the physiologically relevant exposure to hypoxia or by using the chemical DMSO (; ; ; ; ). Because each of these exposures may influence the PHT transcriptome in a manner that is independent of the effect on differentiation, we focused on RNA interactions that are shared by both exposures. We surmised that identification of such interactions may suggest RNA regulatory pathways that govern trophoblast differentiation.
Our informatics-based data offered general insight into changes in the trophoblast transcriptome during differentiation. We noticed that the average number of reads for each lncRNA is about 1.3∼1.5% of the average mRNA read number. Based on Ensembl Gene 87 (), the average length of a lncRNA molecule is 1 kb, and 1.7 kb for mRNA. Therefore, the average expression of lncRNA in PHT cells is between 2.2 and 2.5% of mRNA. This lower expression of lncRNA relative to mRNA was shown in other cell systems (). Note that the reads per library for lncRNAs and mRNAs in the first experimental set are higher because of the larger average library size in that set. Further, our data show that the distribution of the reads per library of “other RNAs” is somewhat different between the two experiments (Figure 2). These two factors likely reflect differences in library size, preparation, and sequencing technology between the two experimental sets. Indeed, even among the “other RNAs,” these factors tend to have a greater effect on smaller RNAs (Supplementary Figure 4). Importantly, these technical factors would not affect our analysis, because we tested the differential expression in the two experimental sets separately, and took library size into account in our statistical tests.
Importantly, we found that lncRNAs, miRNAs, and mRNAs respond to hypoxia and DMSO in different ways. Notably, mRNAs and lncRNAs exhibited a similar change in expression patterns despite of their vastly different expression levels: both RNA types clustered by the experimental conditions (Figure 3), both differentially expressed mRNAs and lncRNAs tended to have larger log2 fold expression change (Figure 4), and both exhibited significant concordance in DE transcripts between the two experimental sets (Figure 5). In contrast, DE miRNA clustering reflected not only the experimental exposure but also the placenta used for PHT cell dispersal and, in the case of hypoxia, even the laboratory performing the sequencing Figure 3). This observation could be explained by the smaller log2 fold expression change in differentially expressed miRNAs (Figure 4), a smaller number of miRNA with a concordant expression change between the two experimental sets (Figure 5), or differences in sequencing technology among the laboratories.
We identified mRNAs that were concordantly differentially expressed across both experimental sets and correlated with altered expression of ≥ 1 lncRNA and miRNA (Supplementary Table 1). Some of the identified transcripts were previously shown to be related to germane trophoblast processes and to complications of pregnancy. SDC1 (Syndecan1) was shown to have a positive correlation with trophoblast differentiation and exhibited lower expression in hypoxia/DMSO exposure compared to standard conditions (). Several identified differentially expressed mRNAs are known to regulate the invasion or adhesion of trophoblasts and other types of cells. KISS1 has a role in early placentation and implantation, ADAM12 (ADAM metallopeptidase domain 12) was shown to control trophoblast fusion through E-cadherin, and NECTIN3 (Nectin Cell Adhesion Molecule 3) is a member of a cell adhesion family of proteins (; ; ). Some of the genes with the most pronounced expression changes are related to preeclampsia, where hypoxia is commonly implicated in disease pathogenesis (). Indeed, reduced placental expression of LGALS13 (Galectin 13), a member of the glycan-binding proteins that regulate innate and adaptive immune responses, was found in women with preeclampsia (). Similarly, the expression of SDC1, NPPB (natriuretic peptide B), GDF15 (growth differentiation factor 15), and ADAMTS6 (ADAM metallopeptidase with thrombospondin type 1 motif 6) all had a lower expression in our experimental exposures, and all are lower in preeclampsia (; ; ; ). In contrast, our data are inconsistent with respect to the expression of PAPPA2, a regulator of trophoblast invasion and migration, or MNDA (myeloid cell nuclear differentiation antigen), which exhibited a lower expression level in our paradigms but is elevated in placentas from women with preeclampsia (; ; ).
A major innovative aspect of our work is the use of the model-based co-expression analysis. The Pearson correlation coefficient is often used as a measurement of co-expressed genes. When the data do not follow multivariate normal distribution, more robust methods, such as the Spearman correlation or the bi-weight mid-correlation, are recommended (). These methods, however, may not be optimal for “real world” data, which are often influenced by confounding variables, such as the sample processing by separate laboratories in our first experimental sets, a factor known to strongly influence expression measurements. Thus, two genes may have significant correlation not because of shared biological pathways, but simply because the expression measurements were performed by a certain laboratory or technology. Using the model-based method proposed in this paper, we can eliminate the “spurious” co-expression caused by known confounding factors. This is critical when data from multiple experiments are combined in a meta co-expression analysis. For example, consider a confounding factor Z with N(0,1) distribution, and two genes X, Y that are conditionally independent, given Z = z, with distribution N(z,1). It then follows that the Pearson correlation between X and Y is 0.5. As mentioned before, under this scenario, the model-based co-expression test is equivalent to testing whether X and Y have a zero partial correlation, given Z. It is easy to see that here, indeed, the partial correlation between X and Y, given Z, is 0. Thus, we correctly capture the independent relation between X and Y given the confounder Z, and eliminate the spurious correlation.
We recognize that this paper has several limitations. First, the data were collected over a period of 4 years, and two laboratories were used to perform RNAseq in the first experimental set. Although we have developed statistical models to address this shortcoming, it might have negatively affected the power of this study. Second, our results were derived from an informatic analysis of the gene expression data on the basis of RNA sequencing libraries. It will therefore require future experimental validation. This is particularly important as the regulation of mRNA by lncRNA and miRNA is a complex process, involving direct and indirect regulation and chromatin remodeling. For the same reason we did not consider the presence of or absence of miRNA-response elements in gene ontology analysis. Finally, our work is limited to in vitro approaches, which may not be fully consistent with changes that occur in vivo. In the future, we plan to collect additional data from our experiments and from public databases to strengthen and validate our results.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA674312; https://www.ncbi.nlm.nih.gov/, PRJNA674329; https://www.ncbi.nlm.nih.gov/, PRJNA674366; https://www.ncbi.nlm.nih.gov/, PRJNA704383; https://www.ncbi.nlm.nih.gov/, PRJNA704399; https://www.ncbi.nlm.nih.gov/, PRJNA704393.
Ethics statement
The studies involving human participants were reviewed and approved by the Institutional Review Board at the University of Pittsburgh. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.
Author contributions
TC, JFM, and YS: study design, manuscript writing. TC and YS: provided study materials. TC and ZC: data generation and analysis. TC, JFM, YO, OB, and YS: data interpretation. All authors contributed to the article and approved the submitted version.
Funding
This project was supported by Eunice Kennedy Shriver National Institute of Child Health and Human Development (NIH/NICHD) grants R01-HD065893 (to YS), Pennsylvania Department of Health Formula Research Funds (to TC and JFM).
Acknowledgments
We thank Elena Sadovsky and Tiffany Coon for technical assistance, Lori Rideout for assistance with manuscript preparation, and Bruce Campbell for editing.
Conflict of interest
YS was a consultant at Illumina, Inc. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.677981/full#supplementary-material
Footnotes
References
1
AghababaeiM.HoggK.PerduS.RobinsonW. P.BeristainA. G. (2015). ADAM12-directed ectodomain shedding of E-cadherin potentiates trophoblast fusion.Cell Death Differ.221970–1984. 10.1038/cdd.2015.44
2
AkenB. L.AchuthanP.AkanniW.AmodeM. R.BernsdorffF.BhaiJ.et al (2017). Ensembl 2017.Nucleic Acids Res.45D635–D642. 10.1093/nar/gkw1104
3
AlsatE.WyploszP.MalassinéA.GuibourdencheJ.PorquetD.NessmannC.et al (1996). Hypoxia impairs cell fusion and differentiation process in human cytotrophoblast, in vitro.J. Cell. Physiol.168346–353.
4
ArnholdtH.MeiselF.FandreyK.LöhrsU. (1991). Proliferation of villous trophoblast of the human placenta in normal and abnormal pregnancies.Virchows Arch. B Cell Pathol. Incl. Mol. Pathol.60365–372. 10.1007/bf02899568
5
BartelD. P. (2009). MicroRNAs: target recognition and regulatory functions.Cell136215–233. 10.1016/j.cell.2009.01.002
6
BeharierO.TyurinV. A.GoffJ. P.Guerrero-SantoroJ.KajiwaraK.ChuT.et al (2020). PLA2G6 guards placental trophoblasts against ferroptotic injury.Proc. Natl. Acad. Sci. U.S.A.11727319–27328. 10.1073/pnas.2009201117
7
BenjaminiY.HochbergY. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing.J. R. Stat. Soc. Series B Stat. Methodol.57289–300.
8
BildiriciI.SchaiffW. T.ChenB.MorizaneM.OhS. Y.O’brienM.et al (2018). PLIN2 is essential for trophoblastic lipid droplet accumulation and cell survival during hypoxia.Endocrinology1593937–3949. 10.1210/en.2018-00752
9
BoydJ. D.HamiltonW. J. (1970). The Human Placenta.Cambridge: Heffer.
10
BurtonG. J.JauniauxE. (2018). Pathophysiology of placental-derived fetal growth restriction.Am. J. Obstet. Gynecol.218S745–S761. 10.1016/j.ajog.2017.11.577
11
BurtonG. J.FowdenA. L.ThornburgK. L. (2016). Placental origins of chronic disease.Physiol. Rev.961509–1565. 10.1152/physrev.00029.2015
12
BurtonG. J.WatsonA. L.HempstockJ.SkepperJ. N.JauniauxE. (2002). Uterine glands provide histiotrophic nutrition for the human fetus during the first trimester of pregnancy.J. Clin. Endocrinol. Metab.872954–2959. 10.1210/jcem.87.6.8563
13
CanfieldJ.ArlierS.MongE. F.LockhartJ.VanwyeJ.Guzeloglu-KayisliO.et al (2019). Decreased LIN28B in preeclampsia impairs human trophoblast differentiation and migration.FASEB J.332759–2769. 10.1096/fj.201801163R
14
Carlevaro-FitaJ.JohnsonR. (2019). Global positioning system: understanding long noncoding RNAs through subcellular localization.Mol. Cell73869–883. 10.1016/j.molcel.2019.02.008
15
CartwrightJ. E.KeoghR. J.Tissot Van PatotM. C. (2007). Hypoxia and placental remodelling.Adv. Exp. Med. Biol.618113–126. 10.1007/978-0-387-75434-5_9
16
ChenB.NelsonD. M.SadovskyY. (2006). N-myc down-regulated gene 1 modulates the response of term human trophoblasts to hypoxic injury.J. Biol. Chem.2812764–2772. 10.1074/jbc.M507330200
17
ChenJ.ThirkillT. L.LohstrohP. N.BielmeierS. R.NarotskyM. G.BestD. S.et al (2004). Bromodichloromethane inhibits human placental trophoblast differentiation.Toxicol. Sci.78166–174. 10.1093/toxsci/kfh046
18
ChenQ.WangY.ZhaoM.HyettJ.Da Silva CostaF.NieG. (2016). Serum levels of GDF15 are reduced in preeclampsia and the reduction is more profound in late-onset than early-onset cases.Cytokine83226–230. 10.1016/j.cyto.2016.05.002
19
ChuT.MouilletJ. F.HoodB. L.ConradsT. P.SadovskyY. (2015). The assembly of miRNA-mRNA-protein regulatory networks using high-throughput expression data.Bioinformatics311780–1787. 10.1093/bioinformatics/btv038
20
DaoudG.SimoneauL.MasseA.RassartE.LafondJ. (2005). Expression of cFABP and PPAR in trophoblast cells: effect of PPAR ligands on linoleic acid uptake and differentiation.Biochim. Biophys. Acta1687181–194. 10.1016/j.bbalip.2004.11.017
21
DeardenL.OcklefordC. (1983). “Structure of human trophoblast: correlation with function,” in Biology of Trophoblast,edsLokeY. W.WhyteA. (Amsterdam: Elsevier).
22
DerrienT.JohnsonR.BussottiG.TanzerA.DjebaliS.TilgnerH.et al (2012). The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression.Genome Res.221775–1789. 10.1101/gr.132159.111
23
DobinA.DavisC. A.SchlesingerF.DrenkowJ.ZaleskiC.JhaS.et al (2013). STAR: ultrafast universal RNA-seq aligner.Bioinformatics2915–21. 10.1093/bioinformatics/bts635
24
DouglasG. C.KingB. F. (1990). Differentiation of human trophoblast cells in vitro as revealed by immunocytochemical staining of desmoplakin and nuclei.J. Cell Sci.96(Pt 1)131–141.
25
DouglasG. C.KingB. F. (1993). Colchicine inhibits human trophoblast differentiation in vitro.Placenta14187–201. 10.1016/s0143-4004(05)80260-7
26
FoxH. (1970). Effect of hypoxia on trophoblast in organ culture. A morphologic and autoradiographic study.Am. J. Obstet. Gynecol.1071058–1064. 10.1016/0002-9378(70)90629-0
27
FrankishA.DiekhansM.FerreiraA. M.JohnsonR.JungreisI.LovelandJ.et al (2019). GENCODE reference annotation for the human and mouse genomes.Nucleic Acids Res.47D766–D773. 10.1093/nar/gky955
28
GandleyR. E.AlthouseA.JeyabalanA.Bregand-WhiteJ. M.McgonigalS.MyerskiA. C.et al (2016). Low soluble Syndecan-1 precedes preeclampsia.PLoS One11:e0157608. 10.1371/journal.pone.0157608
29
GentlemanR. C.CareyV. J.BatesD. M.BolstadB.DettlingM.DudoitS.et al (2004). Bioconductor: open software development for computational biology and bioinformatics.Genome Biol.5:R80. 10.1186/gb-2004-5-10-r80
30
HonC. C.RamilowskiJ. A.HarshbargerJ.BertinN.RackhamO. J.GoughJ.et al (2017). An atlas of human long non-coding RNAs with accurate 5’ ends.Nature543199–204. 10.1038/nature21374
31
HuK. L.ChangH. M.ZhaoH. C.YuY.LiR.QiaoJ. (2019). Potential roles for the kisspeptin/kisspeptin receptor system in implantation and placentation.Hum. Reprod. Update25326–343. 10.1093/humupd/dmy046
32
IyerM. K.NiknafsY. S.MalikR.SinghalU.SahuA.HosonoY.et al (2015). The landscape of long noncoding RNAs in the human transcriptome.Nat. Genet.47199–208. 10.1038/ng.3192
33
JainA.TutejaG. (2021). PlacentaCellEnrich: a tool to characterize gene sets using placenta cell-specific gene enrichment analysis.Placenta103164–171. 10.1016/j.placenta.2020.10.029
34
JauniauxE.PostonL.BurtonG. J. (2006). Placental-related diseases of pregnancy: involvement of oxidative stress and implications in human evolution.Hum. Reprod. Update12747–755. 10.1093/humupd/dml016
35
JiangB.MendelsonC. R. (2005). O2 enhancement of human trophoblast differentiation and hCYP19 (aromatase) gene expression are mediated by proteasomal degradation of USF1 and USF2.Mol. Cell. Biol.258824–8833. 10.1128/mcb.25.20.8824-8833.2005
36
JiangR.WangT.ZhouF.YaoY.HeJ.XuD. (2020). Bioinformatics-based identification of miRNA-, lncRNA-, and mRNA-associated ceRNA networks and potential biomarkers for preeclampsia.Medicine (Baltimore)99:e22985. 10.1097/md.0000000000022985
37
JonesC. J.FoxH. (1991). Ultrastructure of the normal human placenta.Electron Microsc. Rev.4129–178. 10.1016/0892-0354(91)90019-9
38
JunusK.WikströmA. K.LarssonA.OlovssonM. (2014). Placental expression of proBNP/NT-proBNP and plasma levels of NT-proBNP in early- and late-onset preeclampsia.Am. J. Hypertens.271225–1230. 10.1093/ajh/hpu033
39
KlimanH. J.NestlerJ. E.SermasiE.SangerJ. M.StraussJ. F.III (1986). Purification, characterization, and in vitro differentiation of cytotrophoblasts from human term placentae.Endocrinology1181567–1582. 10.1210/endo-118-4-1567
40
KolialexiA.TsangarisG. T.SifakisS.GourgiotisD.KatsafadouA.LykoudiA.et al (2017). Plasma biomarkers for the identification of women at risk for early-onset preeclampsia.Expert Rev. Proteomics14269–276. 10.1080/14789450.2017.1291345
41
KwakY. T.MuralimanoharanS.GogateA. A.MendelsonC. R. (2019). Human trophoblast differentiation is associated with profound gene regulatory and epigenetic changes.Endocrinology1602189–2203. 10.1210/en.2019-00144
42
LangmeadB.TrapnellC.PopM.SalzbergS. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome.Genome Biol.10:R25. 10.1186/gb-2009-10-3-r25
43
LevyR.SmithS. D.ChandlerK.SadovskyY.NelsonD. M. (2000). Apoptosis in human cultured trophoblasts is enhanced by hypoxia and diminished by epidermal growth factor.Am. J. Physiol. Cell Physiol.278C982–C988. 10.1152/ajpcell.2000.278.5.C982
44
LivakK. J.SchmittgenT. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method.Methods25402–408. 10.1006/meth.2001.1262
45
LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.Genome Biol.15:550. 10.1186/s13059-014-0550-8
46
McCarthyC.CotterF. E.McelwaineS.TwomeyA.MooneyE. E.RyanF.et al (2007). Altered gene expression patterns in intrauterine growth restriction: potential role of hypoxia.Am. J. Obstet. Gynecol.19670.E1–70.E6. 10.1016/j.ajog.2006.08.027
47
MouilletJ. F.ChuT.NelsonD. M.MishimaT.SadovskyY. (2010). MiR-205 silences MED1 in hypoxic primary human trophoblasts.FASEB J.242030–2039. 10.1096/fj.09-149724
48
MouilletJ. F.DonkerR. B.MishimaT.CronqvistT.ChuT.SadovskyY. (2013). The unique expression and function of miR-424 in human placental trophoblasts.Biol. Reprod.89:25. 10.1095/biolreprod.113.110049
49
NelsonD. M.JohnsonR. D.SmithS. D.AntebyE. Y.SadovskyY. (1999). Hypoxia limits differentiation and up-regulates expression and activity of prostaglandin H synthase 2 in cultured trophoblast from term human placenta.Am. J. Obstet. Gynecol.180896–902. 10.1016/s0002-9378(99)70661-7
50
NeumanR. I.Alblas Van Der MeerM. M.NieboerD.SalehL.VerdonkK.KalraB.et al (2020). PAPP-A2 and inhibin A as novel predictors for pregnancy complications in women with suspected or confirmed preeclampsia.J. Am. Heart Assoc.9:e018219. 10.1161/jaha.120.018219
51
OhS. Y.ChuT.SadovskyY. (2011). The timing and duration of hypoxia determine gene expression patterns in cultured human trophoblasts.Placenta321004–1009. 10.1016/j.placenta.2011.09.010
52
PalazzoA. F.LeeE. S. (2015). Non-coding RNA: what is functional and what is junk?Front. Genet.6:2. 10.3389/fgene.2015.00002
53
PardiG.CetinI.MarconiA. M.LanfranchiA.BozzettiP.FerrazziE.et al (1993). Diagnostic value of blood sampling in fetuses with growth retardation.N. Engl. J. Med.328692–696. 10.1056/nejm199303113281004
54
PardiG.MarconiA. M.CetinI. (2002). Placental-fetal interrelationship in IUGR fetuses–a review.Placenta23(Suppl. A)S136–S141. 10.1053/plac.2002.0802
55
PrakashG. J.SumanP.GuptaS. K. (2011). Relevance of syndecan-1 in the trophoblastic BeWo cell syncytialization.Am. J. Reprod. Immunol.66385–393. 10.1111/j.1600-0897.2011.01017.x
56
R Development Core Team (2012). R: A Language and Environment for Statistical Computing.Vienna: R Foundation for Statistical Computing.
57
ReymondN.BorgJ. P.LecocqE.AdelaideJ.Campadelli-FiumeG.DubreuilP.et al (2000). Human nectin3/PRR3: a novel member of the PVR/PRR/nectin family that interacts with afadin.Gene255347–355. 10.1016/s0378-1119(00)00316-4
58
RimonE.ChenB.ShanksA. L.NelsonD. M.SadovskyY. (2008). Hypoxia in human trophoblasts stimulates the expression and secretion of connective tissue growth factor.Endocrinology1492952–2958. 10.1210/en.2007-1099
59
RohC. R.BudhrajaV.KimH. S.NelsonD. M.SadovskyY. (2005). Microarray-based identification of differentially expressed genes in hypoxic term human trophoblasts and in placental villi of pregnancies with growth restricted fetuses.Placenta26319–328. 10.1016/j.placenta.2004.06.013
60
SadovskyY.JanssonT. (2015). “Placenta and placental transport function,” in Knobil and Neill’s Physiology of Reproduction,4th Edn, edsPlantT. M.ZeleznikA. J. (San Diego, CA: Academic Press).
61
SahaS.AinR. (2020). MicroRNA regulation of murine trophoblast stem cell self-renewal and differentiation.Life Sci. Alliance3:202000674. 10.26508/lsa.202000674
62
SchaiffW. T.CarlsonM. G.SmithS. D.LevyR.NelsonD. M.SadovskyY. (2000). Peroxisome proliferator-activated receptor-gamma modulates differentiation of human trophoblast in a ligand-specific manner.J. Clin. Endocrinol. Metab.853874–3881. 10.1210/jcem.85.10.6885
63
SchootsM. H.GordijnS. J.ScherjonS. A.Van GoorH.HillebrandsJ. L. (2018). Oxidative stress in placental pathology.Placenta69153–161. 10.1016/j.placenta.2018.03.003
64
ShengF.SunN.JiY.MaY.DingH.ZhangQ.et al (2019). Aberrant expression of imprinted lncRNA MEG8 causes trophoblast dysfunction and abortion.J. Cell. Biochem.12017378–17390. 10.1002/jcb.29002
65
SimonM. C.KeithB. (2008). The role of oxygen availability in embryonic development and stem cell function.Nat. Rev. Mol. Cell Biol.9285–296. 10.1038/nrm2354
66
SoaresM. J.IqbalK.KozaiK. (2017). Hypoxia and placental development.Birth Defects Res.1091309–1329. 10.1002/bdr2.1135
67
StatelloL.GuoC. J.ChenL. L.HuarteM. (2021). Gene regulation by long non-coding RNAs and its biological functions.Nat. Rev. Mol. Cell Biol.2296–118. 10.1038/s41580-020-00315-9
68
TalR. (2012). The role of hypoxia and hypoxia-inducible factor-1alpha in preeclampsia pathogenesis.Biol. Reprod.87:134. 10.1095/biolreprod.112.102723
69
ThanN. G.BaloghA.RomeroR.KárpátiE.ErezO.SzilágyiA.et al (2014). Placental Protein 13 (PP13) - A placental immunoregulatory galectin protecting pregnancy.Front. Immunol.5:348. 10.3389/fimmu.2014.00348
70
ThirkillT. L.DouglasG. C. (1997). Differentiation of human trophoblast cells in vitro is inhibited by dimethylsulfoxide.J. Cell. Biochem.65460–468.
71
Uszczynska-RatajczakB.LagardeJ.FrankishA.GuigóR.JohnsonR. (2018). Towards a complete map of the human long non-coding RNA transcriptome.Nat. Rev. Genet.19535–548. 10.1038/s41576-018-0017-y
72
WagnerP. K.OtomoA.ChristiansJ. K. (2011). Regulation of pregnancy-associated plasma protein A2 (PAPPA2) in a human placental trophoblast cell line (BeWo).Reprod. Biol. Endocrinol.9:48. 10.1186/1477-7827-9-48
73
WakelandA. K.SoncinF.Moretto-ZitaM.ChangC. W.HoriiM.PizzoD.et al (2017). Hypoxia Directs human extravillous trophoblast differentiation in a hypoxia-inducible factor-dependent manner.Am. J. Pathol.187767–780. 10.1016/j.ajpath.2016.11.018
74
XieL.MouilletJ. F.ChuT.ParksW. T.SadovskyE.KnoflerM.et al (2014). C19MC microRNAs regulate the migration of human trophoblasts.Endocrinology1554975–4985. 10.1210/en.2014-1501
75
YusufK.SmithS. D.LevyR.SchaiffW. T.WyattS. M.SadovskyY.et al (2001). Thromboxane A(2) limits differentiation and enhances apoptosis of cultured human trophoblasts.Pediatr. Res.50203–209. 10.1203/00006450-200108000-00007
76
YusufK.SmithS. D.SadovskyY.NelsonD. M. (2002). Trophoblast differentiation modulates the activity of caspases in primary cultures of term human trophoblasts.Pediatr. Res.52411–415. 10.1203/00006450-200209000-00018
77
ZhangB.HorvathS. (2005). A general framework for weighted gene co-expression network analysis.Stat. Appl. Genet. Mol. Biol.4:Article17. 10.2202/1544-6115.1128s
Summary
Keywords
trophoblast, differentiation, hypoxia, lncRNA, miRNA, gene ontology, RNA network
Citation
Chu T, Mouillet J-F, Cao Z, Barak O, Ouyang Y and Sadovsky Y (2021) RNA Network Interactions During Differentiation of Human Trophoblasts. Front. Cell Dev. Biol. 9:677981. doi: 10.3389/fcell.2021.677981
Received
08 March 2021
Accepted
29 April 2021
Published
03 June 2021
Volume
9 - 2021
Edited by
Geetu Tuteja, Iowa State University, United States
Reviewed by
Sha Sun, University of California, Irvine, United States; Michelle Holland, King’s College London, United Kingdom
Updates
Copyright
© 2021 Chu, Mouillet, Cao, Barak, Ouyang and Sadovsky.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Tianjiao Chu, tchu@mwri.magee.eduYoel Sadovsky, ysadovsky@mwri.magee.edu
This article was submitted to Developmental Epigenetics, a section of the journal Frontiers in Cell and Developmental Biology
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.