Transcriptomes flux and change dramatically across both space and time throughout all organismal developmental programs, and during the execution of all normal cellular and organismal functions as well as in disease processes. In human pregnancy, dynamic regulation of the transcriptome occurs in the mother, as well within the fetus (Gao et al., 2018). These trends are reflected in the maternal circulation (Gomez-Lopez et al., 2019) and in the amniotic fluid (Tarca et al., 2020) throughout pregnancy. The transcriptome is a critical effector, and reflection, of the successful initial establishment of the prenatal maternal environment, and serves as both a multidimensional set of markers and a corresponding set of biological effectors for normal biological growth, and parturition. We recently found that prenatal cell-free fetal DNA (cfDNA) is associated with inflammation-induced cytokine release response in a murine model of preterm birth (Gomez-Lopez et al., 2020). After conducting a pilot study (Tarca et al., 2019b) to benchmark the Human Transcriptome Arrays (HTA 2.0) for identification of gestational age- and parturition-related expression changes in the coding and non-coding transcriptome, we have subsequently used this platform to profile a larger patient population at six time points throughout pregnancy in order to establish a baseline reference atlas of gestational age modulated gene expression. Furthermore, microbial genomes in fetal plasma (Witt et al., 2020) were recently detected from umbilical cord blood, highlighting the utility of next-generation sequencing to detect intra-amniotic infection and for personalized diagnostics of high-risk complications of preterm birth. Despite these key advances, little remains known about regulation of the non-coding transcriptome during pregnancy, including small (such as micro) and long non-coding RNAs, and the epigenetic as well as post-transcriptional, cis as well as trans regulatory effects of pregnancy-responsive ncRNAs on the proteins and at the key genomic loci that are functionally relevant to pregnancy and/or its complications. Many non-coding RNAs (lncRNAs, miRNAs, expressed pseudogenes, and other classes of non-coding transcripts) have been shown or implicated to play key roles in active regulatory contexts, albeit–especially for most of the still poorly characterized lncRNAs–to an uncertain extent. Nevertheless, it is by now generally appreciated that the non-coding regions of metazoan genomes exert substantial levels of functional control over cellular activities, do so dynamically under various biological conditions, and are hence hypothesized, and likely, to have direct and specific impacts on the transcriptomic and molecular-pathway responses (Dunham et al., 2012), including from early to term normal pregnancy. Although changes with gestational age of protein-coding genes were confirmed by targeted approaches (Tarca et al., 2019b), in this study, we aimed at characterizing using targeted approaches some of the most highly modulated ncRNAs. The roles of these genes will subsequently need to be functionally validated with direct follow-up gain- and/or loss-of-function studies of these ncRNAs in primary cell culture or organoid models of direct relevance to pregnancy.

RNA from serial maternal blood samples was profiled using microarrays and expression was quantified at whole-transcript level as well as at the individual-exon level. Of the 614 transcript clusters with expression above background and significant correlation with gestational age at blood draw, 104 were initially assigned to the non-coding category (Gomez-Lopez et al., 2019). The longitudinal profiles of these non-coding RNAs belonged to four broad patterns: Up over time, Down over time, Up then Down, and Down then Up. Among the genes most highly modulated during gestation, there were nine long non-coding RNAs (lncRNAs) and long intergenic non-coding RNAs (lincRNAs) (AL355711, BC039551, JHDM1D-AS1, A2M-AS1, NR_034004, LINC00649, LINC00861, LINC01094, MANEA-AS1) (Table 1) and six transcribed pseudogenes (UBBP4, FOXO3B, MKRN9P, LOC441455, PSME2P2, YBX3P1) (Table 2), as well as three classical non-coding regulatory RNAs (MIR4439, SNORD41, and SCARNA2) (Table 3).

We previously reported that transcripts (predominantly, of protein-coding genes) that are differentially expressed with gestational age in normal pregnancy were enriched in genes located on chromosome 14. The differentially expressed non-coding RNA genes, however, were fairly evenly distributed across the human genome (Tables 1, 2, 3). This might imply more diverse lncRNA etiologies of developmental processes that are not tied to a single chromosome-14 immune response locus.

We used the Functional Annotation of the Mammalian Genome (FANTOM)5 Consortium dataset (Forrest et al., 2014) to interrogate the gestational age-associated (TDE) lncRNAs based on their specificity to 1,000 human organs, tissue types, and primary cell cultures in the FANTOM single-base-resolution promoterome resource, accessed through FANTOM5’s ZENBU interface. The majority of these TDE lncRNAs are robustly expressed in the circulating immune cells, as would be expected, given that these RNAs were isolated from maternal peripheral blood. Supplementary File S5 contains all FANTOM5 ZENBU Browser supplementary figure illustrations of the corresponding lncRNA gene structures.

We found that expression profiles of the TDE lncRNAs and lincRNAs modulated during gestation in normal pregnancy are typically tissue-specific, as is expected of functional lncRNAs, and were broadly consistent with pregnancy-associated functions. This supports the contention that these ncRNAs are relevant to the physiology of GA progression in the maternal-fetal environment, rather than being false positives or expression byproducts of GA-responsive protein-coding loci.

Summarily, expression profiles determined by single-stranded 5′Cap Analysis of Gene Expression (CAGE) RNAseq in approximately 1,000 human tissues and cell types represented in the FANTOM5 ZENBU Promoterome Expression dataset (Forrest et al., 2014) indicate strongly specific expression of TDE lncRNAs in immune cells, in particular lymphocytes and monocytes. This is consistent with direct functional roles of these non-coding RNAs in maternal blood, as opposed to purely being biomarkers of GA with functions in other tissues. As transcriptomic biomarkers (Forrest et al., 2014), these lncRNAs, along with their co-expressed mRNAs, define the transcriptome response to pregnancy over advancing GA (Table 1). Additional functional studies are needed in order to determine conclusively whether these TDE lncRNAs have specific regulatory roles in gestation, fetal growth, and/or at the maternal-fetal interface.

Because of the inherent false positive and false negative findings in genome-scale expression microarray experiments, and because of the genomic landscape complexity as well as the comparatively low expression levels of lncRNA genes, a de facto standard emerging in the lncRNA transcriptomics field over the past 2 decades centers on experimental validation of microarray-based lncRNA differential expression results using an independent method in the same samples. Quantitative real-time PCR (RT-qPCR) typically serves as such an independent method. The generally lower expression levels of lncRNAs compared to protein-coding mRNAs make lncRNAs less suitable overall for sybrGreen RT-qPCR. Accordingly, we subjected our top TDE lncRNAs to Taqman RT-qPCR validation, because this method performs better with low-expression targets while the sequence-specificity of the fluorescent Taqman probe guards against false positives. Of the 126 samples selected for RT-qPCR experiments, for 116 there was enough RNA available and analyses were successfully performed.

Of the six lncRNAs tested by RT-qPCR for significant association with gestational age, five were successfully validated (Supplementary Files S3, S4). For four of them (LINC01094, LINC00861, BC039551, and AL355711) RT-qPCR expression was significantly correlated with microarray data, based on which the hypothesis was generated (all p < 0.05) (Figure 2). For LINC01094, LINC00861, BC039551, and A2M−AS1, a direct association with gestational age was also established based on linear mixed effects models in which the gestational age was modeled as a function of gestational age at blood draw and random effects were allowed to account for within patient repeated measurements (p < 0.05) (Figure 3).

Recent advances in sequencing and global unbiased analyses of tissue and timecourse transcriptomic signatures can help reveal the functions of lncRNAs in gene networks, including those in the cell and tissue types most pertinent to pregnancy and maternal-fetal medicine. The high tissue-specificity of TDE lncRNAs in pregnancy would make them potentially attractive therapeutic targets, if the deregulation of any specific TDE lncRNAs were to be discovered during disorders or complications of pregnancy. Expression profiling of these TDE lncRNAs, which we have now characterized in the normal pregnancy GA timecourse, should be also studied in pregnancy complications.

A search of Pubmed-indexed publications for the co-occurrence of terms “gestational age” and “lncRNA” results in a list of 49 publications, as of the time of this writing. However, to date, only four of those publications have investigated the transcriptome (Finn et al., 2014; Gormley et al., 2017; Hussey et al., 2020; Wang et al., 2020), instead of focusing on specific-gene case studies. Specific components of the human lncRNA transcriptome in gestation have recently been elucidated. RNA sequencing to compare normal pregnancy with fetal growth restriction (FGR) (Wang et al., 2020) pinpointed RP11-552M6.1, LINC01291, and ASGR1 (umbilical cord blood) and SFRP2, miR-432-5p, and miR-1306-3p (maternal peripheral blood). Of these non-coding targets functionally associated with and active in FGR (Wang et al., 2020) one was, a Gencode lncRNA originally found, pertinently, in fetal liver in 1999 (Genbank accession number AF090889, found in our UCSC Genome Browser-assisted manual annotation of this locus). RP11-552M6.1 is antisense to the protein-coding gene Kruppel-like factor 12 (KLF12). KLF12 is ubiquitously expressed in numerous organs and tissues, including placenta and uterus (RNA-Seq Data from GTEx, data not shown) and few reported roles in pregnancy. KLF12 is a significant quantitative trait locus from multiple GWAS (Genome-Wide Association Studies) for a host of non-pregnancy associated outcomes such as cancers (Petersen et al., 2010), lipids (Kathiresan et al., 2007), eye disorders (Macgregor et al., 2010), sudden cardiac arrest (Aouizerat et al., 2011), and heart function abnormalities (Sotoodehnia et al., 2010). GWAS is a highly reliable, robust, and well-established method that utilizes population genetics and genomic epidemiology to identify new contributing candidate genes, in an unbiased manner, for common diseases that have a genetic component; these genes can be functionally validated in the laboratory for development into biomarkers or drug targets (Manning et al., 2020). The antisense intronic lncRNA of the KLF12 gene, RP11_552M6.1, is a promising functional target, because KLF12—despite its promiscuous expression—is known to repress AP-2, the alpha-a developmental specific transcription factor. We can envision a role for this locus in FGR of pregnancy that involves tissue- and developmental stage-specific antisense-driven and regulation of the ubiquitously expressed KLF12 gene by the lncRNA. Similarly, we might posit a role for MIR1306 and miR-432-5p in post-transcriptional processing affecting mRNA stability in abnormal pregnancy outcomes, especially for MIR1306, because it is highly expressed in reproductive organs, including the uterus and vagina (data not shown, GTEx).

Despite the growing number of studies of lncRNA in the disorders of pregnancy, the circulating lncRNAome of normal pregnancy over gestational age in humans has been comparatively little investigated. Our study is the first to analyze the lncRNAome over time during normal pregnancy. In addition to GA-dependent lncRNAs in maternal circulation, our analysis, has uncovered other regulatory ncRNAs that were differentially expressed with GA in normal pregnancy in maternal blood: MIR4439 (log ₂ FC =0.41, adj.p= 0.0021, p= 0.0002), which was originally discovered in human B-cells (Jima et al., 2010), SNORD41 (log ₂ FC =0.41, adj. p= 0.0099, p= 0.0011), and SCARNA2 (log ₂ FC = −0.32, adj. p = 0.0044, p = 0.0005). GA was associated with opposite trends in expression levels of the 3 DE regulatory ncRNAs, suggesting selective effects on the network effectors specific for each regulatory ncRNA throughout the timeframe of gestation. For instance, MIR4439 expression dropped in expression roughly mid-gestation and slowly rose again over time until birth. On the other hand, SNORD41 rose subtly but significantly over time, while the lincRNA SCARNA2 (non-coding RNA) dropped rather dramatically over the entire course of pregnancy.

Several of TDE lncRNA genes are robustly expressed in circulating CD14+ monocytes, a cell type strongly functionally associated with preterm labor (Gervasi et al., 2001; Pique-Regi et al., 2019), and other pregnancy complications (Gervasi et al., 2001; Gervasi et al., 2002; Oggé et al., 2010). LINC00649 (chr21:35303516-35343487) (log ₂ FC= −0.35, adj. p = 1.55E-07, p = 3.21E-09) is expressed in immune cell types (as are several other TDE lncRNAs we found), including CD14+ monocytes, CD4+ and CD8+ T-cells, and NK cells. No apparent literature links LINC00649 to any aspects of pregnancy. LINC00649 was shown to be an adverse prognostic marker in AML and is a putative part of a competing endogenous RNA (ceRNA) network (Guo et al., 2020) that includes key AML regulators.

These preliminary results warrant a comprehensive in-silico analysis and experimental validation of miRNA targets, the ncRNA-protein interactions, and their corresponding gene networks associated with pregnancy outcomes. The limited information currently available on LINC00649 complicates predictive analyses (since we do not know which networks encoded at and regulated by this locus are actually utilized during pregnancy). LINC00649 directly harbors a GWAS SNP, rs35184820, significantly associated with menarche (age at onset) in 242,000 European-ancestry individuals (Kichaev et al., 2019). Other GWAS SNPs reside within and near this lncRNA–but not at ATP5O, which alludes to direct functional significance of the lncRNA itself rather than a bystander effect from ATP5O–are reported for hematocrit, erythrocyte counts, glomerular filtration rate, and FEV/FEC ratio. LINC00649 shares a bidirectional promoter with, and (depending on the splice isoform and on the transcription start site used) is antisense to, ATP5O, which encodes an F-type ATPase mitochondrial matrix protein associated with type II diabetes (T2D) in the context of its skeletal muscle mRNA expression profile (Rönn et al., 2009). This gene is also reported among the top 20 genes upregulated in benign epithelial ovarian cysts compared to normal ovarian tissue [see Table 2; (Liu et al., 2017)] but no prior work investigated its role over GA during pregnancy in humans.

LINC00861 (log ₂ FC = −0.40, adj. p = 1.70E-07) is another lincRNA differentially regulated over GA that displays a strong immune cell expression profile, as does LINC00649, including expression in T cells, NK cells, and other immune cell types, according to RNAseq data from the GTEx Portal multiple tissue expression database. Thus, the detection of this lincRNA in maternal blood is expected, and the lincRNA might be relevant to the complications of pregnancy. Intriguingly, the region immediately adjoining this lincRNA transcriptional unit harbors the GWAS SNP rs16901004, significantly associated with endometriosis in 171 European-ancestry cases and 2,934 European-ancestry controls (Sobalska-Kwapis et al., 2017). Previous reports have implicated this lincRNA in survival in ovarian cancer (Zheng et al., 2020), recurrence of hepatocellular carcinomas (Lv et al., 2018), as a clinical marker in breast cancer (Zhang et al., 2017), and in non-smoking COPD (Qian et al., 2018). A role in gestational-age dependent physiology is slightly perplexing in light of the lincRNA’s strong predilection for associating with cancer-related phenotypes and as a biomarker in oncological clinical disease (Zhang et al., 2017; Lv et al., 2018; Zheng et al., 2020). In view of the classical “cancer/stem cell” similarity of the oncotranscriptome with the transcriptome of pluripotent cells and broadly of development, one possibility is a role for the lincRNA in growth and matrix stability in cell and various tissues, but not as a marker for actual carcinogenic and transformative process per se. Additional studies are necessary to clarify the functional basis of our characterization of this lincRNA as a maternal circulating transcriptome biomarker of GA.

The forkhead box FO3B pseudogene (FOXO3B) (log ₂ FC = 0.66, adj. p = 2.24E-15) is a GA-associated transcribed retropseudogene whose genomic context and parental-gene functions are suggestive of a role in pregnancy. FOXO3B is highly expressed in whole blood (like many of our other TDE, GA-responsive lncRNAs and pseudogene transcripts) and in neutrophils (Supp Fig. 11, 25). The FOXO transcription factor genes (FOXOs), which include this pseudogene’s parental gene, are highly associated with pregnancy and its adverse outcomes in humans [for a review see Brosens et al. (2009)], while in mice the FOXOs are required for growth of the preimplantation embryo as a consequence of their involvement in apoptotic processes and cell cycle arrest (Kuscu et al., 2019). The balance of expression levels between the 1 and 3a FOXO transcription factors is protective against oxidative damage during endometrial decidualization (Kajihara et al., 2006). These pivotal roles for FOXOs are expected, given that the conserved “forkhead” family of transcription factors [those containing a forkhead box (FOX) DNA binding domain] includes the FOXOs, and their deregulation originally empowered the discovery of the Drosophila melanogaster mutant phenotype (fkh) that subsequently led to one of the first identifications of an important developmental family of transcription factors [for a review see Hannenhalli and Kaestner (2009)]. The parental gene is DE with GA as well. The FOXO family of transcription factors have long been known as regulators of the cellular changes that accompany maternal tissue growth and remodeling throughout the gestation process (Kajihara et al., 2006; Brosens et al., 2009). Based on its expression pattern, the FOXO3B pseudogene is one of the most striking examples of the DE genes with GA in this study, being present in many maternal reproductive tissues according to GTEx expression profiles- including ovary, fibroblasts, endocervix, uterus, ectocervix, fallopian tubes, EBV-transformed lymphocytes, vagina, testes and mammary tissue- a highly specific expression profile consistent with a function in gestation. We are the first to report that this pseudogene is also differentially expressed with GA during pregnancy. The combination of this pseudogene’s TDE status in pregnancy and the prominence of the female reproductive system for the expression profile strongly alludes to a yet-to-be-characterized function in pregnancy for this transcribed pseudogene.

In addition to the FOXO3B pseudogene, we identified several GA-responsive (TDE) pseudogenes expressed in maternal circulation during normal pregnancy: UBBP4, MKRN9P, LOC441455 (also a makorin pseudogene), PSME2P2, and YBX3P1. The past decade has seen the emergence of multiple experimentally validated modalities of pseudogene-dependent gene regulation, including networks where pseudogene transcripts function as ceRNAs, which are typically non-coding or pseudogene transcripts that compete with their sequence-cognate mRNAs for binding miRNAs and regulatory factors (Salmena et al., 2011; Bosson et al., 2014; Qu et al., 2015), an initially speculative hypothesis which has garnered progressively more extensive support (Vitiello et al., 2020). The net effect of ceRNA networks is, typically, a sequestration of free miRNAs, through miRNA direct binding to cognate binding sites localized in pseudogene transcripts instead of the 3′UTRs of the cognate mRNAs, thereby de-inhibiting the original miRNA-driven post-transcriptional downregulation of those particular mRNAs. The transcription of the corresponding lncRNA or pseudogene is a prerequisite for this process. That we report a large number of TDE, GA-responsive genes in normal pregnancy which fall into this category highlights the heretofore-unrecognized, and potentially extensive, functions of pseudogene-mediated transcriptional networks in the processes of gestation and parturition, and potentially in pregnancy complications in both the mother and fetus.

The lncRNA JHDM1D-AS1, another TDE GA-responsive hit in our analysis, is antisense to the well-characterized protein-coding gene JHDM1D, which encodes a histone demethylase that demethylates H3K9me2 and H3K27me2, and that binds the classic active-promoter, transcription start site histone modification H3K4me3 (Horton et al., 2010; Qi et al., 2010; Tsukada et al., 2010). Reduced expression levels of JHDM1D in the placenta accompany pre-eclampsia (Luo et al., 2018), but a function for this locus during normal pregnancy, if any, remains unknown. RNAi with siJHDM1D, to mimic the reduction in expression levels seen with pre-eclampsia, induced apoptosis in the HTR-8/SVneo trophoblast cell line (Luo et al., 2018), and may increase methylation at HLA-G where JHDM1D binds (Tang et al., 2015; Luo et al., 2018), resulting in a concomitant drop in HLA-G as well, hypothetically also during pre-eclampsia. This pattern correlates with the reduced expression of HLA-G proteins, an immune change found commonly in maternal serum and placental tissues during pre-eclampsia (Yie et al., 2004). We have identified the expression of this antisense lncRNA in amniotic fluid, and we have RT-qPCR-validated the expression of this lncRNA in amniocytes isolated from cultured membranes (A.-S. Goustin, manuscript in preparation), indirectly supporting a likely specific role for this antisense transcript in pregnancy. Here, we report for the first time that JHDM1D-AS1 is dynamically regulated over time during gestation in maternal peripheral blood, slowly rising over increasing GA and slightly prior to birth (log ₂ FC = 0.63, adj. p = 1.62E-09).

LncRNAs are, as a rule, globally expressed at lower levels than are mRNAs of protein-coding genes, while their tissue- and cell-type-specificity equal or exceed those of protein-coding genes and can be used as windows into function (Derrien et al., 2012). In agreement with those fundamental early findings of the ENCODE Consortium, here we found that, while expression profiles of TDE lncRNAs and lincRNAs in FANTOM5 Consortium datasets were low for some though not all of the genes, their expression was overall highly tissue- and organ-specific, often in patterns strongly suggestive of a functional rather than purely a blood-biomarker association with pregnancy. The tissue-specificity of some of the GA responders, especially those with relatively high expression levels in one or relatively few tissues, indicates that—if their functions can be proven—they and their downstream targets may prove potentially attractive for therapeutics. Knowing the expected trajectory of gene expression during normal pregnancy will empower disease diagnostics through personalized transcriptomics; RNAseq of individual-patent samples is affordable, will pinpoint significant deviations from the baseline that we have characterized here, and will define the directions of the desired therapeutic adjustments (up or down, and specific GA) for each gene. Summarily, here we have provided the first ever atlas of lncRNA, lincRNA, antisense lncRNA, and transcribed-pseudogene differential expression in maternal peripheral blood over the timecourse of normal human pregnancy, advancing the field beyond mere transcriptomic biomarkers and toward the discovery of direct pregnancy roles of the lncRNAome. Generating testable predictions is easier for antisense lncRNAs and expressed pseudogenes than for true intergenic lincRNAs, thanks to the known roles of their overlapping sense protein-coding genes and parental protein coding-genes respectively. Therefore, these categories of lncRNAs serve as a valuable starting point in adjusting our focus from a Universe of circulating biomarkers to a constellation of new therapeutic targets for the prevention of pregnancy complications and adverse neonatal outcomes (Huang et al., 2018).