Multimodal immune phenotyping of maternal peripheral blood in normal human pregnancy

Changes in maternal immunity during pregnancy can result in an altered immune state, and as a natural perturbation, this provides an opportunity to understand functional interactions of the immune system in vivo. We report characterization of maternal peripheral immune phenotypes for 33 longitudinally sampled normal pregnancies, using clinical measurements of complete blood counts and major immune cell populations, as well as high parameter flow cytometry for 30 leukocyte antigens characterizing 79 cell populations, and monitoring of 1305 serum proteins using the SomaLogic platform. Cellular analyses characterized transient changes in T cell polarization and more persistent alterations in T and B cell subset frequencies and activation. Serum proteomic analysis identified a potentially novel set of 7 proteins that are predictive of gestational age: DDR1, PLAU, MRC1, ACP5, ROBO2, IGF2R, and GNS. We further show that gestational age can be predicted from the parameters obtained by complete blood count tests and clinical flow cytometry characterizing 5 major immune cell populations. Inferring gestational age from this routine clinical phenotyping data could be useful in resource-limited settings that lack obstetric ultrasound. Overall, both the cellular and proteomic analyses validate previously reported phenotypic immunological changes of pregnancy and uncover potentially new alterations and predictive markers. Research Article Immunology Reproductive biology


Introduction
Epidemiological observations show that maternal immunity is perturbed during normal human pregnancy. Increased risk of severe disease is reported in pregnancy for defined infections, including listeriosis, malaria, and varicella (1)(2)(3). Worse outcomes of influenza infection, and diminished antibody responses to influenza vaccination, are observed in pregnant women (4,5). On the other hand, specific maternal autoimmune disorders, such as rheumatoid arthritis and multiple sclerosis, may regress during pregnancy, whereas symptoms of systemic lupus erythematosus and myasthenia gravis may worsen (6)(7)(8). The immunological changes during pregnancy that underlie these effects are complex and remain poorly understood. Immunological accommodation of the semiallogeneic fetus is mediated primarily by the dedicated tissues of maternal decidua and fetal placenta, which exhibit multiple features of opposing selective forces in viviparity, including imprinting (9). Maternal peripheral immunity is also clearly altered, with changes in systemic IL-4 levels and frequencies of Th2 cells described in normal pregnancy over 25 years ago (10,11). Subsequent work has continued to define the changes in peripheral immune phenotypes to add to our understanding of the physiologically significant modifications of maternal immunity that occur during normal pregnancy (12,13). This understanding could address the clinical susceptibilities or resistance to some pathologies during pregnancy or aid the design of maternal immunization strategies.
Given the broad nature of peripheral maternal immunological changes during pregnancy, high dimensional phenotyping is particularly appropriate. Systems biology analyses of intracellular signaling responses measured by mass cytometry, and combined with other modalities, have been used to define molecular Changes in maternal immunity during pregnancy can result in an altered immune state, and as a natural perturbation, this provides an opportunity to understand functional interactions of the immune system in vivo. We report characterization of maternal peripheral immune phenotypes for 33 longitudinally sampled normal pregnancies, using clinical measurements of complete blood counts and major immune cell populations, as well as high parameter flow cytometry for 30 leukocyte antigens characterizing 79 cell populations, and monitoring of 1305 serum proteins using the SomaLogic platform. Cellular analyses characterized transient changes in T cell polarization and more persistent alterations in T and B cell subset frequencies and activation. Serum proteomic analysis identified a potentially novel set of 7 proteins that are predictive of gestational age: DDR1, PLAU, MRC1, ACP5, ROBO2, IGF2R, and GNS. We further show that gestational age can be predicted from the parameters obtained by complete blood count tests and clinical flow cytometry characterizing 5 major immune cell populations. Inferring gestational age from this routine clinical phenotyping data could be useful in resource-limited settings that lack obstetric ultrasound. Overall, both the cellular and proteomic analyses validate previously reported phenotypic immunological changes of pregnancy and uncover potentially new alterations and predictive markers.
clocks of gestational age (GA) (14,15). In addition to biological understanding of the changes during gestation, these clocks could be used to detect pathologies of pregnancy, such as preterm birth. To complement these data, we have conducted a multimodal high dimensional analysis of longitudinally sampled individuals in the context of normal pregnancies. Here, we identified changes in maternal peripheral immunity, focusing on immune cell populations and the serum proteome. This analysis allowed us to validate previous findings, identify potentially novel changes throughout pregnancy, and determine models predictive of GA that use minimal sets of serum proteins or routinely acquired clinical phenotyping data.

Results
Study design. For 33 pregnancies with vaginal delivery of healthy singletons, 4 peripheral blood samples were obtained from 3 time points throughout gestation as well as once after delivery. The distribution of sampling time points and cohort characteristics are summarized in Figure 1 and Table 1. For all 132 samples immune cell phenotypes were characterized by high parameter flow cytometry, and the serum proteome was analyzed using the SomaLogic platform. Measurements were also obtained from complete blood count (CBC) tests and low parameter clinical flow cytometry quantifying 5 major immune cell populations.
Immune cell population phenotypes. Peripheral blood immune cell phenotypes were studied by high parameter flow cytometry to quantify 79 cellular subsets (Supplemental Tables 1 and 2; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.134838DS1). Populations that significantly changed in frequency were identified by longitudinal paired tests comparing between visit 1 (early gestation), visit 3 (late gestation), or visit 4 (postpartum). We detected 32 population subsets that differed significantly in at least 1 of these comparisons (fold change > 1.2 or < 0.8, q < 0.05) ( Figure 2A). In some cases, the same population was identified by more than one of the parallel staining panels used, and technical validation was provided by detection of similar changes in these populations. An example is the increase in HLA-DR + CD4 + cells after parturition, which was detected by 3 panels as populations 7.1, 7.2, and 34. Among the 32 populations, some changed during gestation and rebounded after parturition, so no difference was observed between early gestation and postpartum, representing transient gestation-associated perturbations. Other populations did differ between the early gestation and postpartum time points, indicating persistent perturbations because of either gestation or parturition ( Figure 2A).
Transient gestation-associated changes in T cell polarization. A transient bias in T cell polarization was observed during gestation that resolved rapidly after parturition. Notably, Th1 and Th17 cell frequencies decreased during gestation, between visits 1 and 3 ( Figure 2A). This occurred not only for CD4 + T helper cells generally but also for CD4 + CXCR5 + T follicular helper (Tfh) populations. The Tfh compartment showed significantly decreased frequencies of cells with type 1/17 phenotype, defined by CXCR3 and CCR6 expression, and significantly increased frequencies of type 2 cells lacking both of these markers ( Figure 2, B and C). The shift in polarization was also observed for CD8 + T cells, with significantly decreased frequencies of CXCR3 + CCR6type 1 cytotoxic T cells during gestation ( Figure 2A). CXCR3 expression by CD8 + T cells can enable recruitment to sites of inflammation, for example by virus-specific cells in acute infection (16). Almost all the T cell polarization phenotypes that changed significantly during gestation, between visits 1 and 3, changed significantly in the opposite direction between visits 3 and 4, which spanned parturition, so significant differences did not remain between the early gestation and postpartum time points (Figure 2A).
Persisting perturbation postparturition. In contrast to the rapidly resolving polarization changes, other differences in T and B cell populations were observed between the extreme time points of our study, the early gestation and postpartum visits 1 and 4 ( Figure 2A). Strikingly, the longitudinal profile of these persistent perturbations differed between T and B cell populations. T cell subsets showed few changes between early and late gestation, other than in polarization, but after parturition numerous populations differed compared with either time point during gestation. During this period, CD8 + T cells skewed from naive CCR7 + C-D45RA + cells to terminal effector CCR7 -CD45RA + or effector memory CCR7 -CD45RAcells ( Figure 2A). Changes in activation differed between CD4 + and CD8 + T cells, with the frequency of HLA-DR + activated CD4 + cells increasing and CD38 + activated CD8 + cells decreasing (Figure 2A). In contrast to these T cell changes, B cell populations that differed between visits 1 and 4 did show significant change before the end of gestation. Frequencies of transitional and activated naive B cells decreased between visits 1 and 3, before rebounding more strongly during parturition, resulting in frequencies significantly higher postpartum compared with early gestation (Figure 2, A and D). Together these observations delineate impacts of pregnancy on the immune system that persist from early gestation to beyond 10 weeks after parturition.
Serum protein characterization. Peripheral blood serum proteins were measured using the SomaLogic platform to quantify 1305 proteins, in samples matching those used for flow cytometry. Longitudinal paired tests comparing between early and late gestation (visits 1 and 3) identified 434 proteins that differed significantly (q < 0.05). We observed no correlations between the cell populations and serum proteins that were identified to change during this period (Supplemental Figure 1). This indicates that the circulating proteins and cell phenotypes behave differently and represent independent data sets. We also detected no evidence that either the cell populations or serum proteins identified to change during gestation differed significantly between individuals in correlation with parity, maternal age, prior miscarriage, duration of gestation, body mass index, or blood pressure (Supplemental Table 3). This is likely due to comparison of a small number of pregnancies with little variation in parameters such as duration of gestation.
Two other studies have recently been reported that used a similar approach for serum proteomic measurement throughout normal pregnancy (17,18). Romero et al. found that 10% of the proteins analyzed changed in abundance as a function of GA, and Aghaeepour et al. used elastic net (EN) modeling to identify 74 proteins that could predict GA, as well as a reduced subset of 8, which was similarly effective (17,18). We first set out to validate these findings, by testing whether the reported protein sets could be used to predict GA, using the measurements for these proteins in our cohort. Excluding 4 proteins that failed quality control parameters in our assay, we used the 70-and 8-protein sets Aghaeepour et al. defined to generate EN models predicting GA. The models were trained on 80% of randomly selected samples and tested on the remaining 20%. For both protein sets, highly significant correlations were seen between observed and predicted GA (r > 0.9, P = 0.0004) (Figure 3, A and B; and Supplemental Figure 2). Thus our results validate the serum proteins reported to predict GA in an independent data set.
Identification of more clinically tractable predictors of GA. Given the predictive analyses possible using serum protein measurements, and the changes in cell population frequencies observed throughout gestation, we sought to predict GA using more accessible features, such as a minimal protein set or routinely available clinical data. Excluding the 70 proteins Aghaeepour et al. defined, we repeated generation of an EN model and found that prediction of GA could still be observed (r > 0.9, P = 0.0004) ( Figure 3C). Seven proteins selected with highest average weight were all significantly different from null models: DDR1, PLAU, MRC1, ACP5, ROBO2, IGF2R, and GNS ( Figure 3D). In pairwise comparisons with the set of 8 proteins Aghaeepour et al. identified, the 7 proteins newly identified showed some correlations but also substantial differences, indicating that these protein sets behave somewhat differently and do not directly mark each other ( Figure 3E). From this comprehensive serum analysis, we therefore identify a novel set of 7 proteins that can be quantified to predict GA.
Given the widespread changes throughout gestation that we observed in cell population frequencies by high parameter flow cytometry, we investigated whether GA could be predicted from immune population frequencies determined at the much lower resolution of standard clinical phenotyping. We applied the same EN modeling approach used for predictive analysis with proteomic data, using all 33 subjects and the 3 time points during gestation. Models were generated using routine clinical phenotyping measurements from CBC tests and clinical flow cytometry that additionally quantified 5 major populations of T, B, and NK cells (Supplemental Table 4). Using 19 parameters from both CBC tests and clinical flow cytometry, models predicting GA were identified with significant correlation between predicted and observed GA (r = 0.45, P = 0.0004) ( Figure 3F). Features that associated with GA were increase in absolute counts of monocytes, decrease in B and NK cell absolute counts, as well as decrease in hemoglobin level ( Figure 3G). We did detect significant changes in total monocyte and NK populations by high parameter flow cytometry, where specifically CD14 + CD16 dim monocytes and activated NK cells were altered most as a fraction of their parent populations (Figure 2A). Changes in peripheral NK cell frequencies throughout gestation are not well established in previous studies, although CD56 bright NK cells show extraordinary enrichment to dominate the decidual leukocyte population at implantation; these cells are transcriptionally distinct from circulating NK cells (19)(20)(21). An increase with GA has previously been reported for the total monocyte population (20). This population can be quantified by a CBC test without clinical flow cytometry, and indeed, when using only 9 parameters measured from the CBC test alone, models predicting GA could still be identified where monocyte counts and hemoglobin levels were the dominant features (r = 0.35, P = 0.0012) (Supplemental Figure 3). These results support the idea that simple, cost-effective, routine clinical phenotyping measurements could be used to predict GA.

Discussion
We report high dimensional analyses characterizing immune cell populations and the serum proteome, in maternal peripheral blood sampled longitudinally during normal pregnancy. In both the immune populations and serum proteome, we validate previously reported phenotypic changes, and add new details and markers to the changes in maternal peripheral immunity that have been identified. Expanding the changes defined adds to our understanding of the physiologically significant modifications of maternal immunity that occur during normal pregnancy and may aid the design of maternal immunization strategies.
During gestation, transient T cell polarization from Th1/17 to Th2 was observed, not only for CD4 + T helper cells generally but also for CXCR5 + Tfh cells specifically and also CD8 + T cell subsets. T cell polarization has long been described to alter during pregnancy, with changes documented in Th1/2/17 subsets of the broad T helper peripheral CD4 + population, although sometimes with contradictory results (10,22,23). Our data support a peripheral polarization bias from Th1/17 to Th2 but that the effects are transient, with the postpartum time point no different from early gestation. Our novel observations that Tfh and even CD8 + cells show this polarization could potentially modify how these polarization shifts are thought to alter responses to viral infections and autoantigens, which have often focused on Th cells (24)(25)(26). Tfh cells are important for vaccine responses, and the changes in this population may be particularly relevant for optimization of maternal vaccination (27).
Other cellular phenotypic changes showed differences in comparisons between early gestation and postpartum time points. The dynamics of T cell changes, observed particularly in frequencies of CD8 + cell subsets, and more widely in terms of activation status, was consistent with previous reports that broad CD4 + or CD8 + T cell populations show persistent alteration relative to time points from early gestation, in measurements of cell counts and in vitro cytokine production (20). Although immune activation has been implicated in the onset of parturition, the observation that more than 10 weeks after parturition CD4 + T cell activation is still increased compared with our early gestation time point could potentially affect responses in diverse scenarios, such as infection, autoimmunity, and pathologies associated with chronic inflammation (28)(29)(30). It is possible the changes in cell populations could affect transfer of passive immunity to breastfed neonates because CD8 + T cells particularly have been observed to survive in intestinal Peyer's patches of nursed infants (31). B cell populations showed a different profile of changes, with population frequencies and activation states decreasing during gestation before rebounding more strongly after parturition, to result in levels that were higher postpartum compared with early gestation. Our observations are consistent with previous reports that transitional B cells increase within days after parturition (32). Changes throughout the B cell compartment could derive from endocrine factors because hormones such as human chorionic gonadotropin have been implicated in modifying B cell differentiation and function (33). These changes in B cells could contribute to mechanistic understanding of some of the specific alterations in immunity that are observed during pregnancy, such as reduced response to influenza vaccination (4). In addition to changes in population frequencies, it will be important to understand alterations in cellular functional responsiveness during pregnancy because intracellular signaling markers at endogenous levels and in response to stimulation have been found to dominate modeling of GA (14). We applied an EN modeling approach to validate sets of 70 and 8 serum protein sets previously reported to predict GA and to identify a novel set of 7 serum proteins that can be used to predict GA. The serum protein marker of GA with highest average weight was epithelial discoidin domain-containing receptor 1 (DDR1), which is known to be involved in mammary gland development (34). DDR1 has previously been associated with birthweight in the context of disparate genome-wide approaches. Analyses of selection in primate evolution, and differential methylation between pre-and full-term birth in cord blood, both implicated the region encoding DDR1 (35,36). Others among our identified markers of GA have also been previously associated with pregnancy. Urokinase-type plasminogen activator (PLAU) is involved in extravillous trophoblast migration during placentation, detected to be upregulated in normal pregnancy, and associated with pregnancy disorders (37)(38)(39). Cation-independent mannose-6-phosphate receptor, also known as insulin like growth factor 2 receptor (IGF2R), was one of the first genes identified to be imprinted in mice, although in humans only a minority of individuals show biased expression from the maternal allele (40,41). This modeling of GA is restricted to the period from 11 to 37 weeks for which samples were collected; it would be of significant interest to characterize earlier time points of gestation with prospectively collected cohorts.  A and B) EN models generated from both a set of 70 proteins, and a subset of 8 proteins previously identified to predict GA, demonstrated significant correlation between observed and predicted GA when using serum proteomic data from 3 time points during pregnancy for the 33 women in our study. Predicted GA is shown with mean and standard deviation of 100 model iterations randomly sampling training and test sets, with red dashed line indicating a linear regression with 95% confidence interval. (C) Excluding these 70 proteins and selecting from the remaining 1124 serum proteins measured, EN models still showed significant correlation between observed and predicted GA. (D) For the 7 proteins newly identified to predict GA, frequency of selection and average weights differed significantly from null models and are plotted in red showing mean and standard deviation from the 100 model iterations randomly sampling training and test sets. In gray are shown null models with mean and standard deviation from random permutation of sample labels 25 times and 100 random samplings of the training and test sets. (E) The relationship between new and previously identified proteins predicting GA was assessed by pairwise Pearson's correlations of protein relative intensities during gestation in our data set. Eight proteins previously identified to predict GA by Aghaeepour et al., and the 7 proteins identified in this study, are numbered 1-8 and 9-15, respectively, and defined as follows: 1, chorionic somatomammotropin hormone; 2, prolactin; 3, macrophage colony-stimulating factor 1 receptor; 4, polymeric immunoglobulin receptor; 5, proto-oncogene tyrosine-protein kinase receptor Ret; 6, glypican-3; 7, granulins; 8, α-fetoprotein; 9, macrophage mannose receptor 1; 10, tartrate-resistant acid phosphatase type 5; 11, N-acetylglucosamine-6-sulfatase; 12, cation-independent mannose-6-phosphate receptor; 13, epithelial discoidin domain-containing receptor 1; 14, urokinase-type plasminogen activator; and 15, roundabout homolog 2. (F) Standard clinical phenotyping measurements from CBC tests and clinical flow cytometry can also be used to predict GA. For the 33 subjects sampled at 3 time points during pregnancy, 19 routinely acquired clinical parameters were used to generate EN models for which predicted and observed GA correlated significantly, with black dashed line marking equal GA. (G) Frequency of selection and average weights are shown for the significant model features, plotted as in D.
In addition to the set of 7 serum proteins, we also show that GA can be predicted from the routinely acquired clinical data of CBC tests and clinical flow cytometry for 5 major immune populations. The major features selected in these models were monocyte counts and hemoglobin levels. These models may be investigated in pregnancies with pathologies, to potentially shed light on the mechanisms conferring risk. Preterm labor, for example, has been shown to be predicted by cell-free RNA, white blood cell counts, and corticotrophin-releasing hormone, but premature rupture of membranes remains the most accurate indicator of delivery within 48 hours (42,43). Although much larger studies will have to be done before these observations can be adapted clinically, the ability to infer GA from standard CBC tests and clinical flow cytometry measurements could be most valuable in resource-poor settings. In these circumstances obstetric ultrasound is often not possible to monitor fetal growth for detection of abnormalities that indicate increased risk of neonatal mortality (44).

Methods
Subjects and clinical phenotyping. Subjects were included in this study with a pregnancy of GA between 10 and 20 weeks, estimated by ultrasound or last menstrual period, and maternal age between 18 and 45 years. Exclusion criteria were a hemoglobin reading of less than 8 g/dL or any medical condition adversely affecting the immune system or requiring immunomodulating medications. In total 33 individuals were studied, of which 76% described their ethnicity as white and 21% reported prior miscarriage. Subjects were analyzed at 4 time points, with clinical phenotyping, high parameter flow cytometry, and serum proteomic analysis performed for all 132 samples. Clinical phenotyping was performed in the NIH clinical center and comprised CBC tests and clinical flow cytometry, together measuring counts and frequencies for platelets, CD4 + T cells, CD8 + T cells, CD19 + cells, NK cells, neutrophils, and monocytes, as well as hemoglobin levels (Supplemental Table 4).
High parameter flow cytometry. PBMCs were isolated by separation with Ficoll-Paque PLUS and cryopreserved according to Center for Human Immunology (CHI) protocols (45). High parameter flow cytometry was performed using the Human Immunology Project Consortium (HIPC) panels as previously described (46,47). Briefly, 5 parallel 10-color panels with a total of 30 unique markers enable detection of 79 subsets of PBMCs represented as a fraction of their parent population (Supplemental Tables 1 and 2). Staining was performed using the HIPC lyophilized antibody plates, each with up to 10 samples in addition to controls, after an additional incubation with LIVE/DEAD Fixable Blue Dead Cell Stain (Thermo Fisher Scientific). Acquisition was performed with a Becton Dickinson LSRFortessa, using DIVA 8 software, acquiring 250,000 cells for each sample. Subsequent analysis to determine population frequencies used FlowJo version 9.6.2. Compensation performed with unstained cells and Becton Dickinson compensation beads was used to aid acquisition monitoring. Subsequently a final compensation matrix was calculated using FlowJo during postacquisition analysis. Gates used to define all 79 populations are shown for a representative sample (Supplemental Figures 4-8). For 1 sample, subject 008 at visit 4, high parameter flow cytometry data were not obtained.
Serum proteomic analysis. Peripheral blood serum was isolated using SST tubes and cryopreserved according to CHI protocols (45). Serum proteomic analysis used the SOMAscan 1.3k Assay (SomaLogic). This is an aptamer-based assay able to detect 1305 protein analytes, optimized for analysis of human serum (48,49). Briefly, aptamers are short single-stranded DNA sequences modified to confer specific binding to target proteins and can be highly multiplexed for discovery of biomarker signatures. The proteins quantified include cytokines, hormones, growth factors, receptors, kinases, proteases, protease inhibitors, and structural proteins. A complete list of analytes measured can be found at http://somalogic.com/wp-content/uploads/2017/06/SSM-045-Rev-2-SOMAscan-Assay-1.3k-Content.pdf. The assay was performed according to manufacturer specifications, with data then inspected using a web tool and subjected to quality control procedures as previously described, after which results from 1194 somamers were used for analyses (50,51).
Data and code availability. The data used in this study, comprising CBC tests, clinical flow cytometry, high parameter flow cytometry, and serum proteomic analysis by SomaLogic, as well as the R code and documentation to reproduce the results of our analyses, are available at https://github.com/kotliary/pregnancy.
Statistics. Statistical analysis was performed using R/Bioconductor. Longitudinal changes in PBMC populations or serum proteins were evaluated by Wilcoxon's signed-rank paired test, and P values were corrected for multiple comparisons with Benjamini-Hochberg method. Modeling of GA using serum proteins or cell populations was performed by EN method using the eNetXplorer package varying parameter α between 0 (ridge regression) and 1 (lasso regression) with step 0.1 (52). For each α the regularization parameter λ was selected by a crossvalidation procedure optimizing correlation between observed and predicted GA. For this crossvalidation the modeling was run 100 times, each time using a random 80% of samples as the training set and the remaining 20% as the test set. The best parameter α was chosen for models with highest average correlation between observed and predicted GA. To estimate significance of model coefficients, eNetXplorer computes null models, randomly permuting sample labels 25 times, each time with 100 random samples of training and test sets as for the model.
Study approval. For this study volunteers were enrolled in 10-I-0205, an NIH protocol approved and monitored by the National Institute of Allergy and Infectious Diseases/NIH Institutional Review Board (Bethesda, Maryland, USA). All subjects provided informed consent prior to their participation in the study.

Author contributions
LB, SMO, HZ, and RS recruited subjects and collected samples. JC, A Biancotto, and A Babyak performed experimental assays. RA, YK, FC, and KLH analyzed the data. RA, YK, YB, and JST wrote the manuscript. SMH, YB, JST, and CSZ conceived and oversaw the study.