The Cellular Transcriptome in the Maternal Circulation During Normal Pregnancy: A Longitudinal Study

Pregnancy represents a unique immunological state in which the mother adapts to tolerate the semi-allogenic conceptus; yet, the cellular dynamics in the maternal circulation are poorly understood. Using exon-level expression profiling of up to six longitudinal whole blood samples from 49 pregnant women, we undertook a systems biology analysis of the cellular transcriptome dynamics and its correlation with the plasma proteome. We found that: (1) chromosome 14 was the most enriched in transcripts differentially expressed throughout normal pregnancy; (2) the strongest expression changes followed three distinct longitudinal patterns, with genes related to host immune response (e.g., MMP8, DEFA1B, DEFA4, and LTF) showing a steady increase in expression from 10 to 40 weeks of gestation; (3) multiple biological processes and pathways related to immunity and inflammation were modulated during gestation; (4) genes changing with gestation were among those specific to T cells, B cells, CD71+ erythroid cells, natural killer cells, and endothelial cells, as defined based on the GNF Gene Expression Atlas; (5) the average expression of mRNA signatures of T cells, B cells, and erythroid cells followed unique patterns during gestation; (6) the correlation between mRNA and protein abundance was higher for mRNAs that were differentially expressed throughout gestation than for those that were not, and significant mRNA-protein correlations were observed for genes part of the T-cell signature. In summary, unique changes in immune-related genes were discovered by longitudinally assessing the cellular transcriptome in the maternal circulation throughout normal pregnancy, and positive correlations were noted between the cellular transcriptome and plasma proteome for specific genes/proteins. These findings provide insights into the immunobiology of normal pregnancy.


INTRODUCTION
Pregnancy represents a unique immunological state in which the immune system of the mother undergoes adaptations that allow her to tolerate the semi-allogenic conceptus (1)(2)(3). Indeed, pregnancy is divided into three different immunological stages based on cytokine profiles (4). Pioneer studies indicated that, while the innate immune system is upregulated to protect the mother against infection and the fetus from rejection (5,6), the adaptive immune response toward paternal/fetal antigens seems to be selectively suppressed [i.e., driven toward a T-helper (Th)2like phenotype] (7)(8)(9)(10)(11). Specifically, the cellular components of the innate immune system in the maternal systemic circulation are activated as evidenced by increased numbers of monocytes and granulocytes (12,13). Such innate immune cells display an activated phenotype, comparable to that observed in women with sepsis (13), and exhibit enhanced functionality (phagocytosis, respiratory burst activity, and cytokine production) (14)(15)(16)(17)(18). The humoral components of the innate immune system are also boosted during pregnancy (5). For example, complement components and acute phase proteins are increased in the circulation of pregnant women (19)(20)(21)(22)(23)(24). In contrast to the innate immune system, the cellular (e.g., T cells and B cells) and humoral (e.g., antibodies) components of the adaptive immune system in the maternal circulation during normal pregnancy have received less attention.
The systemic intravascular inflammatory response during normal pregnancy is especially activated in women who experience the physiological process of labor at term (18,25) and in those who undergo pregnancy complications such as preterm labor (18,26,27), preterm premature rupture of membranes (28), and preeclampsia (13,17,(29)(30)(31)(32). Therefore, the systemic immune response reflects both physiological and pathological processes, and the early detection of these changes may lead to the discovery of non-invasive biomarkers for obstetrical disease.
Herein, for the first time, we aimed to provide a roadmap of the modulations in the cellular transcriptome in maternal circulation during normal pregnancy. In addition, we assessed whether gestational mRNA changes of the cellular transcriptome correlate to those of the plasma proteome during normal pregnancy.

Study Design
We conducted a prospective longitudinal study that enrolled women attending the Center for Advanced Obstetrical Care and Research of the Perinatology Research Branch, NICHD/NIH/DHHS; the Detroit Medical Center, and Wayne State University School of Medicine. Based on this cohort, we designed a retrospective study that included 49 women with normal pregnancy who delivered at term and had 4-6 blood samples collected throughout gestation [median number of samples = 5, interquartile range (IQR) = 5-6] (n = 282). Blood samples were collected at the time of a prenatal visit, scheduled at 4-week intervals from the first or early second trimester until delivery in the following gestational age intervals: 8-<16, 16-<24, 24-<28, 28-<32, 32-<37, and >37 weeks. All patients provided written informed consent and the use of biological specimens, as well as clinical and ultrasound data, for research purposes were approved by the Institutional Review Boards of Wayne State University and NICHD. All experiments were performed in accordance with relevant guidelines and regulations.

RNA Extraction
RNA was isolated from PAXgene R Blood RNA collection tubes (BD Biosciences, San Jose, CA; Catalog #762165), as described in the PAXgene R Blood miRNA Kit Handbook. Purified RNA was quantified by UV spectrophotometry using the DropSense96 R Microplate Spectrophotometer (Trinean, Gentbrugge, Belgium), and quality was assessed by microfluidics using the RNA ScreenTape on the Agilent 2200 TapeStation (Agilent Technologies, Wilmington, DE, USA).

Microarray Analysis
RNA was processed and hybridized to GeneChip TM Human Transcriptome Arrays 2.0 (P/N 902162) according to the Affymetrix GeneChip TM WT Pico Reagent Kit Users Guide (P/N 703262 Rev. 1) as follows: Biotinylated cDNA were prepared from 20-50 ng total RNA. Labeled cDNA were hybridized to the arrays in a GeneChip TM Hybridization Oven 640 by rotating at 60 rpm, 45 • C for 16 h. Arrays were then washed and stained in the Affymetrix Fluidics Station 450 and scanned using the Affymetrix 3000 7G GeneChip TM Scanner with Autoloader. Raw intensity data were generated from array images using the Affymetrix AGCC software.

Data Analysis
Preprocessing Affymetrix Human Transcriptome Arrays CEL files were preprocessed using Robust Multi-array Average (RMA) (33) implemented in the oligo package (34) and annotation from the hta20sttranscriptcluster.db package of Bioconductor (35). Since samples were profiled in several batches as a part of a larger study, correction for potential batch effects was performed using the removeBatchEffect function of the limma package in Bioconductor. After batch effect correction, data from the sample collected at the time of labor from the 21 women who had spontaneous term labor were removed to avoid confounding gestational age-related changes with those due to the onset of labor at term. The final analysis set of 261 transcriptomes was used in downstream analyses described below.

Expression Calling
Transcript clusters (typically one or two per unique gene) were deemed present in a given sample if one of its probesets (targeting a specific exon) was expressed above background (p-value for expression above background p DBAG <0.05) determined using the Transcriptome Analysis Console (version 4.0) (ThermoFisher Scientific). Genes were retained if deemed present in >25% of the 261 samples.

Differential Expression
Expression profiles were visually inspected to determine consistency of the data in sequential samples collected from Frontiers in Immunology | www.frontiersin.org the same woman. One of 261 samples consistently had the lowest value for a large fraction of the genes and was deemed as outlier and removed from further analysis. Linear mixedeffects models (36) were then used to fit log 2 gene expression data as a function of gestational age (continuous) and included cubic polynomial terms of gestational age as fixed effects and a random intercept term for each woman. Significance p-values for the association of gene expression and gestational age were determined using likelihood ratio tests between a model with and without gestational age terms. A False Discovery Rate adjusted p-value (q-value) <0.1 and a fold change (FC) of >1.25 were required for significance. Fold change was determined as the ratio of the highest vs. lowest average expression from 10 to 40 weeks of gestation. Linear mixed-effects models were fit using the lmer function, while the likelihood ratio tests were performed using the anova function available in the lme4 R package (36).

Gene Ontology and Pathway Analysis
Gene ontology and pathway analysis was conducted using a hypergeometric test on Gene Ontology (GO) (37) and Developmental FunctionaL Annotation at Tufts (DFLAT) databases (38), as well as on Curated Gene Sets (C2) collection from the Molecular Signatures Database (MSigDB) database (39). In addition, enrichment tests were performed for tissue specificity and chromosomal locations of genes. Tissue-specific genes were defined as those with median expression 30 times higher in a given tissue than the median expression of all other tissues documented in the Gene Atlas (40) as previously described (41).
Unless otherwise stated, all enrichment analyses were based on a hypergeometric test and accounted for multiple testing with q < 0.05 being considered a significant result. In all enrichment analyses, the background gene list was defined as the compendium of genes deemed present in >25% of the samples. In this analysis, we tested whether previously reported cell-type specific mRNA signatures derived by single-cell RNA-Seq studies of placenta tissues (42) were modulated with advancing gestation in normal pregnancy. The 13 cell types identified by Tsang et al. (42) were: B cells, T cells, monocytes, cytotrophoblasts, syncytiotrophoblast, decidual cells, dendritic cells, endothelial cells, erythrocytes, Hofbauer cells, stromal cells, vascular smooth muscle cell, and extravillous trophoblasts. The mRNA signatures for these cell types were first quantified in each patient sample by averaging expression data over genes part of each signature. Before averaging, the data for each gene was first standardized by subtracting the mean and dividing by standard deviation of expression across term samples (>37 weeks). Cell-type specific expression averages were then fit as a function of gestational age using linear mixed-effects models, as described above for the analysis of data of individual genes.

Assessment of mRNA Protein Correlations
Maternal plasma abundance of 1,125 proteins in 71 samples collected from 16 of the women included in the current study were obtained from the S1 File of Erez et al. (43). The correlation between each mRNA and corresponding protein pair was assessed by fitting linear mixed-effects models with the response being the protein abundance and the predictor being the mRNA expression. These models included a random intercept term to account for the repeated observations from the same subject. The meaning of the mRNA coefficient in this model is change in log 2 protein abundance for one unit change in log 2 gene expression. The significance of the protein-mRNA correlation was assessed by the t-score for the regression line slope, and false discovery rate adjustment of resulting p-values was performed across all mRNAprotein pairs that were tested. A q-value <0.1 was considered a significant result.

Longitudinal Patterns of the Cellular Transcriptome Throughout Normal Pregnancy
The mRNA profiles of longitudinal maternal blood samples were determined at exon level resolution by microarrays. The characteristics of the study population are shown in Table 1.
A total of 26,458 protein-coding mRNA transcript clusters were expressed above background levels in at least 25% of the samples, as were 5,706 non-coding RNA transcript clusters.  being annotated to immune processes. Chromosome 14 includes genes of critical importance for immunity (44); therefore, these data show that pregnancy has a strong effect on the maternal immune system. To define clusters of expression trajectories during gestation, we focused on 112 of the 614 significant transcript clusters that changed more than 1.5-fold from 10 to 40 weeks of gestation. Three distinct clusters of expression modulation emerged: genes that (1) steadily increased throughout gestation (89 genes; Figure 2, red cluster), (2) steadily decreased throughout gestation (12 genes; Figure 2, green cluster), and (3) decreased prior to mid-gestation followed by an increase (11 genes; Figure 2, blue cluster). These results indicate that the expression of the most highly modulated genes increases with advancing gestational age.
Of note, the 19 mRNA transcript clusters (corresponding to 16 unique genes) that changed more than 2-fold in expression during pregnancy all increased from 10 to 40 weeks of gestation (Figure 3, gray lines correspond to individual pregnancies and blue lines show the average expression). The expression of these 16 genes increased from early to late pregnancy and tended to plateau near term, with the exception of 2 genes (interferon-induced protein 27 and 44-like) (Figure 3). Several of these most highly modulated genes are related to host immune response (e.g., MMP8, DEFA1B, and DEFA4) (45,46), again emphasizing the immune response adaptations during normal pregnancy.

Biological Processes, Pathways, and Immune Cell Signatures Associated With Advancing Gestation in Normal Pregnancy
We performed gene ontology enrichment analysis to interpret the changes in gene expression occurring throughout gestation. We identified 157 biological processes modulated during gestation, which included cellular and humoral immunity, defense response, response to external biotic stimulus (e.g., bacteria and viruses), regulation of lipid storage, interleukin-1beta production and secretion, and erythrocyte development, among others ( Table 2). An additional 134 biological processes altered during gestation were found when querying the Developmental FunctionaL Annotation at Tufts (DFLAT) database, such as stress response, FIGURE 2 | Clustering of average gene expression profiles throughout normal pregnancy. Average profiles of genes that change throughout normal pregnancy and have a fold change >1.5 were clustered using hierarchical clustering. The distance metric used in the clustering was 1-Pearson correlation. Three clusters were identified: Cluster 1 (red, 89 genes), Cluster 2 (green, 12 genes), and Cluster 3 (blue, 11 genes). Note that, in this figure, the average gene expression profiles vs. gestational age were reset so that their value is 0 at 10 weeks of gestation.
immune system development, cytokine response, and regulation of angiogenesis ( Table 3).
Enrichment analyses were then expanded to canonical pathways and gene sets from the Molecular Signatures Database (MSigDB), and 53 such pathways were found to be associated with advancing gestation. These included the Reactome database (47) pathways: immune system, adaptive immune system, cytokine signaling in immune system, and immunoregulatory interactions between a lymphoid cell and a non-lymphoid cell, as well as the KEGG database (48) pathways: natural killer cell-mediated cytotoxicity, antigen processing and presentation, and graft vs. host disease ( Table 4).
We then aimed at determining the origin of observed transcriptional activity by using the GNF Gene Expression Atlas to define genes predominantly expressed in specific human tissue or cell types, as previously described (41). This analysis revealed that gene sets specific to CD4+ and CD8+ T cells, CD71+ erythroid cells, CD105+ endothelial cells, and CD56+ NK cells, among others, were over-represented among the mRNAs that were modulated during gestation (q < 0.05) ( Table 5). In addition, genes reported to be specific to fetal organs (liver, lung, and brain) and the placenta were also enriched among significant genes (q < 0.05) ( Table 5). These data show that maternal and fetal cell-specific transcripts found in the maternal circulation are being modulated with advancing gestation.
The average abundance of cell type-specific gene sets recently defined using single-cell transcriptomics (42) were also analyzed for systematic changes with gestational age at blood draw in our cohort. This analysis revealed that expression of mRNAs specific to three cell subtypes were dynamically altered throughout normal pregnancy: (1) the T-cell-specific mRNA signature decreased from the first to second trimester, followed by a subsequent increase during the third trimester ( Figure 4A); (2) the B-cell-specific mRNA signature decreased steadily throughout gestation ( Figure 4B); and (3) the expression of genes specific to nucleated erythroid cells (HBZ, ALAS2, and AHSP) significantly increased as gestation progressed ( Figure 4C). These findings demonstrate, for the first time, that single-cell RNA-Seqderived signatures of erythroid cells change throughout normal gestation in maternal whole blood, while trends found for T cells and B cells were similar to those reported in whole blood (49) and by using cell-free RNA analysis (42).

Correlation Between the Cellular Transcriptome and Plasma Proteome Throughout Normal Pregnancy
Transcription does not always correlate with protein translation (50,51). Therefore, we investigated the correlation between the mRNAs that were modulated throughout gestation and their corresponding protein abundance. Maternal plasma abundance of 1,125 proteins in 71 samples collected from 16 of the women included in the current study was previously reported (43,52). First, we assessed the mRNA-protein correlation for 53 of the 614 transcript clusters that changed throughout gestation and for which abundance data for the corresponding protein were available. These correlations were compared to those of 1,011 mRNA-protein pairs that did not change with gestation. The mRNA-protein correlations were significantly higher for    transcripts that changed throughout gestation compared to those that did not (Wilcoxon test for comparing t-scores of the linear regression slope obtained by linear mixed-effects models for each mRNA-protein pair, p = 0.01) (Figure 5). Among the 53 transcripts that changed throughout gestation, BPI, IGHG1, CXCL10, GNLY, and GZMA had a significant mRNA-protein correlation as assessed by both linear mixed-effects models and a naïve Spearman correlation test (q < 0.05 for both analyses) (Figure 6). Notably, two of these genes (GNLY and GZMA) were included in the T-cell-specific mRNA signature that was modulated overall throughout gestation.

DISCUSSION
Principal Findings of the Study The correlation between mRNA and protein abundance was higher for mRNAs that were differentially expressed throughout gestation than for those that were not (p = 0.01). (8) Significant and positive mRNA-protein correlations (q < 0.05) were observed for BPI, IGHG1, CXCL10, and two members of the T-cell mRNA signature (GNLY, GZMA). The expression trends and variability in expression of individual genes and meta-genes in normal pregnancy (nomograms) derived herein will be the basis for future studies aiming at developing biomarkers for obstetrical disease.

Transcriptomic Changes During Pregnancy
Previous studies have investigated the cellular (53, 54) and cellfree (55) transcriptome in the maternal circulation at different time points during normal pregnancy using either 3-prime-end biased microarrays or targeted approaches. The current study, however, is the first to quantify at exon-level resolution the cellular transcriptome during normal pregnancy in up to six samples per pregnancy. More than one-half (54%, 277/514) of the unique differentially expressed genes identified herein were also among the 2,321 genes (q < 0.1) reported by Heng et al. (54) to change from 17-23 to 27-33 weeks of gestation. Similarly, 47% (242/514) of the genes found in this study were among the 3,830 genes reported by Al-Garawi et al. (53) as changing from 10-18 to 30-38 weeks. The overlap between the genes reported as differentially expressed in these two studies and those identified herein is significant (Fisher's exact test p < 0.0001 for both). However, unlike in the two previous studies involving a pair of samples from each woman, the availability of four to six longitudinal samples per patient in this study enabled us to capture more complex expression trajectories in the window of 10-40 weeks of gestation, and to identify distinct clusters of such gene expression trajectories. Compared to another recent study by Ngo et al. (55) that involved more frequent sampling than used herein, our study has the advantage of an unbiased assessment of the whole cellular transcriptome as opposed to a targeted assessment of genes that are placenta-, immune-, and fetal liver-specific. Of note, among the 14 immune-specific cell-free mRNAs selected by Ngo et al. (55) as best predictors of gestational age at blood draw, 11 were also identified in our study, with CEACAM8, DEFA4, LTF, and MMP8 being among those with highest fold change. Although our results are somewhat consistent with those reported by Ngo et al. (55), it should be noted that cellular and cell-free transcripts can follow   different patterns in similar physiological and pathological processes (56).

Correlations Between the Cellular Transcriptome and the Plasma Proteome Throughout Normal Pregnancy
The finding that the maternal transcriptome features inflammation-related processes and pathways that are being activated in preparation for labor at term is in agreement with our previous studies in gestational tissues [cervix (57), myometrium (58), membranes (59)] and a similar longitudinal study of the maternal plasma proteome (52). In addition to finding several common biological processes that are modulated in both the maternal plasma proteome and cellular transcriptome (such as defense response, defense response to bacterium, defense response to fungus, regulation of bone resorption, leukocyte migration) we assessed, for the first time, the extent of the agreement in whole blood mRNA and protein changes with gestation in the same set of samples. Although the correlations between mRNA and protein expression reported in the literature are notoriously poor, recent studies showed that mRNA-protein correlation is higher for mRNAs that are differentially expressed in a given condition than for those that are not (51). Our finding that the correlation of mRNAprotein pairs is higher for transcripts changing with gestation Frontiers in Immunology | www.frontiersin.org than those who do not is therefore consistent with previous observations (51).

T Cells in the Maternal Circulation During Normal Pregnancy
Maternal T cells are implicated in the physiological processes occurring throughout gestation (60)(61)(62)(63). Effector and activated T cells are found at the maternal-fetal interface before (64)(65)(66)(67)(68)(69)(70) and during (71)(72)(73) spontaneous labor at term, and these cells are associated with the timing of term parturition (74). Effector T cells are also found in the maternal blood prior to Shah et al. (75) and during (76) labor at term. In the current study, we demonstrated that the T-cell-specific mRNA expression in the maternal circulation was decreased prior to mid-gestation but upregulated from mid-gestation until term. Moreover, for two of 19 genes of this signature, there was a significant correlation between cellular mRNA and plasma proteomic profiles; this is consistent with recent cytomic and proteomic studies in the maternal circulation (77,78). In addition, we recently showed the same u-shaped pattern of expression for the T-cell mRNA signature during gestation in a smaller set of patients profiled with RNA-Seq and qRT-PCR platforms (49). Together, these findings illustrate the importance of maternal T cells during normal pregnancy.
Alteration of systemic T-cell populations has also been implicated in preterm parturition (79)(80)(81), especially since aberrant activation of these cells can induce the onset of preterm labor (82,83). On the other hand, the absence of T cells in a mouse model caused an increased susceptibility to endotoxininduced preterm birth, which was reversed by adoptive transfer of CD4+ T cells (84). From a histopathological standpoint, T cells are detected in placental lesions related to maternal anti-fetal rejection such as villitis of unknown etiology (85)(86)(87), chronic chorioamnionitis (88), and chronic deciduitis (89), which have also been linked to the onset of term and preterm labor (88,(90)(91)(92)(93)(94). The chronic nature (95,96) of these lesions has led our group to propose them as indicators of maternal anti-fetal rejection, which can lead to preterm labor or even fetal death (86,90,93,94,(97)(98)(99)(100). Future studies are needed to establish whether the early detection of T-cell alterations in the maternal circulation may identify pregnancies at risk for obstetrical disease such as preterm labor/birth and fetal death.

B Cells in the Maternal Circulation During Normal Pregnancy
Several studies have suggested a role for B cells in the maintenance and success of pregnancy (101)(102)(103)(104)(105)(106). Circulating CD5+ (B1) B cells were shown to decrease during pregnancy, only returning to normal levels after parturition (107). This finding was later shown to occur in mice, where a decreased influx of newly generated B cells to the blood and spleen was observed while mature B cells were increased in uterinedraining lymph nodes (108). An expansion of marginal zone B cells also ensued (108,109), which was proposed to participate in the production of protective antibodies during pregnancy (109,110). Accordingly, maternal serum antibody concentrations increased concomitantly with B-cell population changes (109), possibly as a result of the anti-inflammatory microenvironment maintained at the maternal-fetal interface throughout most of pregnancy (111). Such antibody production has been considered the primary contribution of B cells to maternal-fetal tolerance during pregnancy (101).
Interleukin-10-producing regulatory B cells (Bregs) have also been described as important players during pregnancy (112,113). Such adaptive immune cells increased in normal pregnancy in an hCG-dependent manner (113,114) and suppressed effector T-cell cytokine production (113). Trophoblast cells facilitated the conversion of IL10-deficient B cells into IL10-expressing B cells (114), which is in line with a previous report showing that the adoptive transfer of Bregs restored maternal-fetal tolerance (112). Indeed, pregnant women treated with the B-cell-depleting treatment rituximab had a higher occurrence of pregnancy loss (115), although further investigation of this phenomenon is warranted (116).
In the current study, we showed that the B-cell-specific mRNA signature moderately decreased throughout pregnancy. Our observations correspond to a previous report indicating that the majority of maternal peripheral B-cell subsets are reduced in late gestation compared to the non-pregnant state (117), whereas Bregs are upregulated (117). The combined effects of such dynamic changes on the overall circulating B-cell mRNA signature are therefore minimal, as we have demonstrated here. Taken together, these findings suggest that, while total maternal peripheral B cells are mostly maintained, subset-specific changes occur throughout pregnancy.
Fetal nucleated erythroid cells have been detected in the maternal blood (121,155,156) where they may be a source of cell-free fetal DNA (155). Previous reports showed that neonatal FIGURE 5 | Distribution of mRNA-protein correlation t-scores. The correlation between mRNA and protein abundance was assessed by linear mixed-effects models using data collected from 71 samples provided by 16 women. The distribution of t-scores for the linear correlation slope is shown for 51 mRNAs differentially expressed with gestation and 1011 mRNAs that were not differentially expressed. CD71+ erythroid cells have immunomodulatory functions on cord blood leukocytes (157)(158)(159), and that their direct contact with maternal peripheral immune cells increases the secretion of pro-inflammatory mediators by such cells (159). Therefore, it is likely that the trafficking of CD71+ erythroid cells from the fetus into the mother directly affects maternal immune responses (159). Nucleated erythroid cells have also been described in the placenta, where their presence is correlated with the number of such cells in the cord blood (160), and these cells also display immunomodulatory properties in vitro (161). FIGURE 6 | Correlation between cellular transcriptome and maternal plasma proteome throughout normal pregnancy. Aptamer-based protein abundance measurements are plotted against mRNA expression. R: naïve Spearman correlation coefficient; p: likelihood ratio test p-value from linear mixed-effects models assessing the linear correlation accounting for repeated measurements from the same subjects.
Herein, we showed that the erythroid cell-specific mRNA signature was upregulated throughout gestation in the maternal circulation. This finding is in line with previous reports showing that fetal microchimerism increases during pregnancy (162,163). In addition, a recent study showed that CD71+ erythroid cells are increased in the maternal circulation throughout gestation, peaking during the third trimester and falling to baseline levels after delivery (164). Together, these findings illustrate that erythroid cells, most likely of fetal origin, are present in the maternal circulation and their transcriptome is modulated as gestation progresses. These data provide a possible mechanism whereby the developing fetus can modulate maternal immunity.

CONCLUSION
We have reported herein a detailed characterization of the longitudinal maternal whole blood transcriptomic changes in normal pregnancy. We have shown that these changes are genome-wide, yet we found that chromosome 14 was particularly enriched in genes modulated with advancing gestation. There was a significant overlap in expression changes described herein with those previously reported in whole blood analyses based on only two time points, while some of the most strongly modulated mRNAs identified herein were also previously reported as the best predictors of gestational age in cell-free RNA analyses of maternal blood. Our systems biology approach to the interpretation of these expression changes in the maternal cellular transcriptome during pregnancy revealed significant longitudinal patterns of expression for immune-related gene sets, such as those specific to T cells, B cells, and erythroid cells. Moreover, for the first time, we demonstrated positive correlations between the cellular transcriptome and plasma proteome for specific genes, including those expressed by T cells. The expression trajectories of protein coding and non-coding transcripts in normal pregnancy described herein may serve as references and hence enable the discovery of biomarkers for obstetrical disease.

DATA AVAILABILITY STATEMENT
The raw and summarized microarrays gene expression data are available as a Gene Expression Omnibus series (https://www. ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121974).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Boards of Wayne State University and NICHD. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
AT, RR, SH, and NG-L conceived the research. SH and RR supervised the enrollment of the patients and collection of samples. AT and NG-L carried out research and drafted the manuscript. GB contributed to data visualization and to the preparation of data submission to the Gene Expression Omnibus. AT analyzed data. AT, RR, and NG-L interpreted the data. SH, RR, PP, JK, and SB provided feedback on the manuscript. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We thank Dan Lot and Dr. Susan Land for conducting the RNA extraction at the Applied Genomics Technology Center of Wayne State University in Detroit, Michigan. We acknowledge Dr. Christopher Krebs for conducting the microarray experiments at the DNA Sequencing Core at the University of Michigan in Ann Arbor, Michigan.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.

2019.02863/full#supplementary-material
Supplementary File 1 | Differential gene expression with gestational age. Each row corresponds to one of the 614 transcript clusters associated with gestational age. ID, Affymetrix transcript cluster identifier; SYMBOL, gene symbol; ENTREZ, Entrez gene database identifier; Name, Gene name; mRNA, mRNA identifier; Source, mRNA identifier source; locus.type, coding vs. non-protein coding assignment; Chromosome, chromosome number; strand, chromosome strand; start/stop, genomic coordinate start and stop for the transcript cluster; p-value, significance value for the polynomial relation between log gene expression and gestational age; FC, log 2 fold change representing the log 2 ratio of the highest and lowest fitted value from 10 to 40 weeks of gestation (see Figure 3). The sign is positive if value at 40 weeks is higher than the value at 10 weeks, and negative otherwise. adj.P.Value, False discovery adjusted p-value.
Supplementary Figure 1 | Longitudinal gene expression profiles of genes associated with gestational age. Each figure shows data for one of the 614 transcript clusters associated with gestational age. The y-axis represents the log 2 normalized gene expression, while the x-axis represents gestational age (weeks). Each line corresponds to one patient and each dot to one sample. The thick blue line represents the linear mixed effects model fit. The title in each plot, represents the transcript cluster identifier, gene symbol, gene name, p-value, and fold change. The meanings of these annotations are as in the legend of Supplementary File 1.