Longitudinal profiling of human blood transcriptome in healthy and lupus pregnancy

Systemic lupus erythematosus (SLE), an autoimmune disease characterized by excess type I IFN, predominantly affects women and often presents during their childbearing years. SLE patients are at increased risk of pregnancy complications, including preeclampsia (PE), fetal or neonatal death, and fetal morbidity from growth restriction and preterm delivery (Clowse et al., 2008; Bundhun et al., 2017). These adverse outcomes affect over one fifth of SLE pregnancies (Buyon et al., 2015). Progress was made in characterizing clinical factors and circulating antiangiogenic factors as predictors of poor pregnancy outcomes, but identification of patients destined for complications remains challenging (Buyon et al., 2015).

Pregnancy requires tolerance of fetal alloantigens encoded by paternal genes. A variety of factors at the maternal–fetal interface protect the fetus from maternal immune responses, including expression of HLA-G (Hunt and Geraghty, 2005), expansion and local accumulation of maternal T regulatory cells (Jiang et al., 2014), immunosuppressive trophoblast-derived exosomes, and epigenetic modifications of decidual stromal cells that limit effector T cell trafficking (PrabhuDas et al., 2015). Absence of these regulatory mechanisms leads to inflammation and trophoblast dysfunction with ensuing placental insufficiency, poor fetal growth, preterm birth, and in some cases, the maternal syndrome of PE.

Experimental models of pregnancy in lupus-like conditions implicate inflammation at the maternal–fetal interface as an essential factor in fetal death and growth restriction (Girardi et al., 2003; Berman et al., 2005). Activation of complement, recruitment of neutrophils, and production of proinflammatory cytokines and antiangiogenic factors are drivers of poor placental development in mice treated with antiphospholipid antibodies and in nonautoimmune models of PE (Girardi et al., 2003, 2006; Gelber et al., 2015), and there is evidence that these mechanisms also contribute to adverse outcomes in human pregnancies (Shamonki et al., 2007; Cohen et al., 2011; Kim et al., 2016). In addition, increased neutrophil numbers and neutrophil extracellular traps (NETs) infiltrating placental intervillous spaces associate with inflammatory, vascular changes in patients with SLE and PE (Marder et al., 2016). Neutrophils are a source of interferogenic nucleic acids (Caielli et al., 2016; Lood et al., 2016) and IFNα has been shown in vitro to prevent stabilization of endothelial cells and thus potentially contribute to PE (Andrade et al., 2015). Furthermore, IFNβ lowers the threshold for inflammatory-driven cytokine release systemically and at reproductive sites (Cappelletti et al., 2017).

Transcriptional biomarkers have been used to study placental tissue, uterine endothelial cells, and the endometrium during healthy pregnancy (Gack et al., 2005; Schmidt et al., 2005; Kitaya et al., 2007; Yockey and Iwasaki, 2018). Measuring these biomarkers, however, requires invasive procedures with associated risks. Blood transcriptional profiling has been widely reported in SLE (Baechler et al., 2003; Bennett et al., 2003; Chaussabel et al., 2008; Banchereau et al., 2016) but never before in SLE pregnancy.

Blood profiling during late stages of pregnancy identified PE-associated cell proliferation, differentiation, and apoptosis-related transcriptional changes (Textoris et al., 2013). We hypothesized that a similar approach applied from early pregnancy might identify SLE patients at high risk for complications. Accordingly, we evaluated the blood transcriptome in individuals from the PROMISSE (Predictors of pRegnancy Outcome: bioMarkers In antiphospholipid antibody Syndrome and Systemic Lupus Erythematosus) study. PROMISSE is the largest multicenter, multiethnic, and multiracial study to prospectively assess clinical and laboratory predictors of adverse outcomes in SLE patients with inactive or mild/moderate disease activity at conception. We profiled the blood transcriptome of SLE patients and healthy pregnant women by microarray and conducted multicolor flow cytometry in a subset of them. To identify the earliest molecular changes associated with pregnancy onset, the blood transcriptome of healthy women undergoing assisted reproductive technologies (ARTs) was profiled before and after embryo implantation by RNA sequencing.

Both in healthy and noncomplicated SLE pregnancies, as well as during the peri-implantation period, we observed a striking transcriptional modulation of immune networks associated with SLE pathogenesis. Incomplete down-regulation of several of these pathways was associated with SLE pregnancy complications, supporting their potentially pathogenic role at the maternal–fetal interface.

Longitudinal blood transcriptional profiles were applied to samples from the prospective, multicenter, observational PROMISSE study. Among PROMISSE patients, adverse outcomes occurred in 19.0% of pregnancies, including fetal (4%) or neonatal (1%) death, preterm delivery <36 wk because of placental insufficiency (9%), and small-for-gestational-age neonate (10%; Buyon et al., 2015). Our study included 92 SLE pregnant, 43 healthy pregnant (H-P), 20 SLE nonpregnant (SLE-NP), and 34 healthy nonpregnant (H-NP) women. For pregnant subjects, the average gestational age at recruitment was 12.2 weeks. Samples drawn at specific time points (P1: <16 wk of gestation [WG]; P2: 16–23 WG; P3: 24–31 WG; P4: 32–40 WG; and between 8 and 20 wk postpartum [PP]) were used for microarray analysis (Fig. 1 A). SLE-pregnant patients were further stratified according to pregnancy outcomes: no complications (NC; n = 46), PE (n = 24), and other fetal-related complications in the absence of PE (OC; n = 22). There were no differences in hydroxychloroquine use among patient groups. For patients receiving steroids, the mean daily dose of prednisone was <10 mg. Unequal sample size distribution of patients in P3 and P4 was due to complication-related premature delivery in the SLE-OC and PE groups (Fig. 1 A).

In addition, we studied 25 subjects undergoing ART by blood RNA sequencing (RNA-seq; Fig. 1 B). Successful pregnancy rate for this cohort was 80% (20/25). There were no significant demographic differences between those who conceived compared with those who did not (Table S1). Blood was drawn on the day of embryo transfer, in the midluteal phase (4–5 d after embryo transfer), on the day of pregnancy test (9 d after embryo transfer), and during the first trimester of pregnancy (P0: 4–6 wk after embryo transfer), when applicable.

We first characterized the blood transcriptional changes elicited over time during healthy pregnancy. To identify differentially expressed transcripts (DETs) at various time points, disease status, and/or complications while accounting for individual variation and missing data, we developed a linear mixed model that included both fixed and random effects. Using H-NP women as baseline controls, 9,576 DETs were identified. Hierarchical clustering revealed prominent transcriptional down-regulation, in line with the highly dynamic immunomodulatory features associated with pregnancy (Mor et al., 2017; Fig. 2 A). Principal component analysis highlighted a gradual shift in signatures away from H-NP from <16 WG that stabilized by 16–23 WG through delivery. The healthy PP signature only partially returned to H-NP baseline, indicating that sustained blood transcriptional changes remain detectable several weeks after delivery (Fig. 2 B).

To interpret these signatures, we leveraged a framework of 260 coexpression modules specifically designed for whole-blood transcriptome analysis, which can detect alterations in both leukocyte frequency and population-intrinsic transcriptional signaling (Chaussabel and Baldwin, 2014). To simplify data interpretation, we depict modules M1–7 (97 out of 260), which display the most conserved coexpression behavior across independent datasets. The global module fingerprint of pregnancy, which combines profiles from four time points compared with H-NP, revealed up-regulation of neutrophil (M5.15), myeloid inflammation (M3.2, M4.2), and erythropoiesis (M2.3, M3.1, M6.18) signatures (Fig. 2, C and D). Conversely, antiviral and lymphocyte-related pathways were downregulated, including modules ascribed to IFN responses (M1.2, M3.4, M5.12), B cells (M4.10), plasma cells (M4.11, M7.7, M7.32), T cells (M4.1, M6.15), and natural killer (NK) cell/cytotoxicity (M3.6, M4.15).

To analyze the evolution of these signatures at each quarter, we used Q-Gen (Turner et al., 2015), a gene-set enrichment approach leveraging QuSAGE (Yaari et al., 2013) that integrates statistics from linear mixed models (Table S2). Significantly altered modules (false discovery rate [FDR]–adjusted P ≤ 0.05, absolute expression >30% in any comparison) were represented as a heatmap (Fig. 3 A), with the five most significant also represented as a line graph (Fig. 3 B). Increase in erythropoiesis and inflammation was detectable as early as P1, while enrichment in neutrophils appeared at P2. Modulation of these signatures was sustained until delivery and subsequently diminished at PP to levels lower than those of H-NP controls. Down-regulation of IFN response- and lymphoid lineage–related modules was detected as early as P1. To validate these changes during healthy pregnancy, whole-blood FACS analysis was conducted on 11 healthy pregnant individuals, three of whom were also transcriptionally profiled (Table S3). An increase in absolute numbers of circulating neutrophils (granular CD14⁻ CD66b⁺) and a decrease in transitional B cells (CD20⁺ IgD⁺ CD24⁺ CD38²⁺) during pregnancy (P1–4) were detected, suggesting that part of the observed transcriptional changes reflect fluctuations in the frequencies of specific leukocyte subsets (Fig. 3 C and Table S4).

To capture early molecular changes associated with pregnancy onset, we monitored the blood transcriptome during embryonic implantation in an independent cohort of 25 age-matched healthy women undergoing ART (Fig. 1 B). Using a mixed model adjusting for subject and time point, combined with Q-Gen analysis incorporating blood modules, early systemic transcriptional changes could be detected at pregnancy onset. Modules that significantly changed at one or more time points compared with their own baseline (embryo transfer) in failed or successful pregnancy are represented as a heatmap (Fig. 4 A). The IFN response increased during the midluteal phase (midluteal vs. embryo transfer; fold change [FC] = 2.16; nominal P = 0.03) but was rapidly down-regulated by pregnancy test in successful pregnancy (pregnancy test vs. embryo transfer; FC = 1.67; nominal P = 0.04), while in failed implantations it was further up-regulated at the time of pregnancy test (pregnancy test vs. embryo transfer; FC = 3.14; nominal P = 0.02). Conversely, erythropoiesis (M2.3) was down-regulated at midluteal (midluteal vs. embryo transfer; FC = 0.79; nominal P = 0.05) but subsequently up-regulated at pregnancy test in successful pregnancies (P0 vs. midluteal; FC = 1.5; nominal P = 0.02), in contrast to its down-regulation in failed pregnancies (pregnancy test vs. midluteal; FC = 0.62; nominal P = 0.04; Fig. 4 B and Table S5). Many IFN-regulated genes were up-regulated in failed pregnancies at the time of pregnancy test (nominal P ≤ 0.05, FC >1.25; Fig. 4 C). These results demonstrate that some of the significant systemic transcriptional changes in maternal hematopoietic compartments start as early as embryonic implantation. They highlight the down-regulation of major immune pathways that are linked to SLE pathogenesis, including the IFN response, and suggest that erythropoiesis-related pathways may be linked to pregnancy-induced immune regulation.

To assess the transcriptome of noncomplicated SLE pregnancy (SLE-NC), we compared the signatures of SLE-NP and SLE-NC women with their respective healthy controls using module maps (Fig. 5 A, top and middle). As reported before for adult and pediatric patients (Chiche et al., 2014; Banchereau et al., 2016), SLE-NP displayed several modular differences compared with healthy controls, including increased IFN and neutrophil signatures (M1.2, M3.4, M5.12, and M5.15; Fig. 5 A, top). A plasma cell signature (M4.11) was only mildly increased, likely reflecting the low disease activity in enrolled patients (mean SLE Pregnancy Disease Activity Index [SLEPDAI] score at screening: 2.5; Table S6). The gene expression signature of SLE-NC compared with H-P controls retained elements of the prototype SLE-NP signature, such as IFN and plasma cell–related transcripts, but included additional up-regulated modules such as cell cycle (M3.3; Fig. 5 A, middle). Myeloid inflammation signatures (M3.2, M4.2, M4.13, M4.14, and M7.27), which were differentially expressed in SLE-NP compared with H-NP, were similarly increased in SLE-NC and H-P (Fig. 2 C). Compared with SLE-NP controls, SLE-NC exhibited remarkable conservation of healthy pregnancy patterns, including down-regulation of all three IFN as well as lymphoid and plasma cell modules and up-regulation of neutrophil and erythropoiesis modules (Fig. 5 A, bottom; and Fig. 2 C). Thus, the transcriptome dynamics of uncomplicated SLE pregnancy mirrors that of healthy pregnancy, but several disease-associated signatures remain dysregulated compared with healthy pregnant controls.

We next compared the profiles of SLE-NC and H-P at each time point and observed 3,138 DETs between the two groups at one or more time points (Fig. 5 B). The largest differences were observed at P2 (1,491 transcripts) and PP (1,252 transcripts). Q-Gen analysis of SLE-NC using the SLE-NP baseline as comparator at each time point confirmed that most elements of the H-P signature were reproduced over time in SLE-NC (Fig. 5 C). Despite widespread transcriptional down-regulation, the increase in neutrophil, erythropoiesis, and myeloid inflammation modules observed in H-P between P1 and P2 was conserved in SLE-NC. Normalization of the SLE-NP and SLE-NC profiles to those of H-NP confirmed that the most significant modular patterns of H-P were reproduced in SLE-NC (Fig. 5 D). Interestingly, while the IFN response down-regulation in SLE-NC did not reach H-NP or H-P levels, the plasma cell signature steadily decreased below the SLE-NP baseline from P1, reaching levels comparable to H-P by P4 and returning to SLE-NP levels after delivery. These results support that the modulatory effects on major transcriptional networks that take place during healthy pregnancy are conserved in SLE pregnancy. The steady and consistent decrease in plasma cell transcripts below H-NP down to H-P levels suggests that modulation of this pathway is a key element of successful SLE pregnancy.

Up to 19% of SLE pregnancies in the PROMISSE study resulted in either maternal and/or fetal–perinatal complications. This substudy was enriched in the number of complications by design and included 50% of SLE-PE/OC. Because the pathogenic mechanisms leading to these complications are likely to be different, we classified patients into PE and “complications other than PE” (OC) groups. H-P and SLE-NC women delivered on average at term (39.3 and 38.7 WG, respectively), while delivery time for complicated pregnancies was 30.8 WG (SLE-PE) and 31.5 WG (SLE-OC; Fig. 6 A and Table S6).

To identify maternal blood profiles associated with complications, we compared signatures of SLE-PE or SLE-OC with those of SLE-NC at each time point during pregnancy using a linear model combined with Q-Gen analysis. Significantly altered modules are depicted as a heatmap in Fig. 6 B, and the most differentially expressed modules among SLE groups are represented as line graphs normalized to H-NP in Fig. 6 C. Higher fold expression of the neutrophil signature was detected in pregnancies destined for PE compared with all other groups at P1 (PE vs. NC; FC = 1.44; FDR-adjusted P < 0.0001), which supports data from experimental models (Gelber et al., 2015). T cell (M4.1) and B cell (M4.10) signatures were down-regulated in early stages in all pregnancy groups, but more significantly in PE (Fig. 6 C and Fig. S1). In addition, suppression of IFN, cell cycle, and especially plasma cell responses were less pronounced in SLE-OC compared with PE. The dynamic pattern of each differentially expressed module by patient group is represented in Fig. S1. Interestingly, the kinetics of inflammation, myeloid lineage, and especially erythropoiesis transcriptional pathways were conserved in all SLE groups and H-P. These signatures are therefore mostly affected by pregnancy, regardless of health status.

High SLE disease activity at conception is associated with increased risk of pregnancy complications (Clowse et al., 2005; Buyon et al., 2015). In the PROMISSE cohort, disease activity score measured by SLEPDAI was slightly higher at the first time point (P1) in SLE-PE (4.43) compared with SLE-OC (2.86) and SLE-NC (2.35) (P = 0.026; Table S6). A second measure of disease activity, the Physician Global Assessment, was comparable across all groups of patients (Table S6). To confirm that PE signatures were not due to overrepresentation of subjects with high disease activity, we compared PE and OC with NC by a mixed model using the earliest time point (P1) with or without the initial SLEPDAI score as a covariate. Modular maps (Fig. S2) representing Q-Gen FC in gene expression in PE vs. NC and OC vs. NC displayed more significant differences after adjustment for SLEPDAI score. This analysis specifically highlights the increased neutrophil signature (M5.15) in SLE-PE and the increased plasma cell signature (M4.11, M7.7, M7.32) in SLE-OC. In addition, it highlights an up-regulation of erythropoiesis and protein synthesis/degradation pathways (M2.3, M4.5, M6.16) in SLE-OC. These data support that signatures related to SLE-PE/OC were not confounded by higher disease activity. The SLEPDAI score was measured at baseline and each trimester. Except for baseline, however, it was not obtained simultaneously with samples for transcriptional studies, and therefore it could not be included as a covariate in analyses of later time points.

To complement transcriptomics results, multicolor FACS analysis was performed in a subset of patients (27 SLE-NC, 5 SLE-PE, and 4 SLE-OC) at P1 and P2, when pregnancy complications were not yet clinically apparent. Transcriptional profiles were simultaneously obtained in 10 (28%) of these patients. Among the 84 cell populations/phenotypes assessed, 42 showed significantly altered frequency in any comparison of PE vs. NC, OC vs. NC, and PE vs. OC in a mixed-model analysis (Fig. 6 D and Table S3). PE showed increased frequency of immature neutrophils at this time point (granular fraction, CD14⁻CD16⁻CD11b⁻ cells; nominal P = 0.017; Fig. 6 E). OC patients displayed increased frequency of activated T cells, as defined by inducible costimulatory molecule (ICOS) expression, both in the total CD4⁺ T cell compartment (nominal P = 0.038) as well as within the follicular helper T (Tfh) cell subset (CD4⁺ CXCR5⁺ ICOS⁺; nominal P = 0.0034). Within Tfh cells, we detected an expansion of CXCR3⁺ cells (Tfh1; nominal P = 0.0012). Within the B cell compartment, transitional B cells (CD20⁺ CD24⁺ CD38²⁺) were significantly expanded in OC (nominal P = 0.02). Early expansion of immature neutrophils in PE and activated T cells in OC might underlie the overexpression of neutrophil and T cell signatures observed in these patient groups, respectively (Fig. 6 B). These results suggest that both innate and adaptive arms of immunity are modulated at early stages of SLE pregnancy and that their dysregulation associates with the development of complications. To facilitate access to data and results from our analyses, we generated and formatted results from statistical analysis that can be uploaded into the genBart shiny app in R. In this interactive interface, users can upload our genBart file (downloadable at http://dx.doi.org/10.17632/hg596dwdr6.1) to view sample information, generate transcriptional heatmaps and module maps, access the mixed model and Q-Gen results of the microarray, browse gene lists driving modular expression, access FACS analysis results, and generate plots for subsets of interest. The interface features are detailed in Fig. S3.

Finally, we sought to identify early transcriptional biomarkers that could predict SLE pregnancy–related PE. Machine learning approaches were applied to P1 transcriptional profiles to predict pregnancy outcomes (PE = 21; NC = 45; total = 66) and compared with previously reported laboratory and clinical variables, including lupus anticoagulant, antihypertensive medication use, Physician Global Assessment score >1, low platelet count, race, and antiangiogenic factors (sFlt1, PlGF, and sFlt1/PlGF; Buyon et al., 2015; Kim et al., 2016). A higher percentage of patients with lupus anticoagulant was present in complicated vs. noncomplicated SLE pregnancies, but, likely due to sample size, this variable did not reach significance (25% in PE, 38% in OC, and 9% in NC; Fisher’s exact test, P = 0.08; Table S6).

We trained the model using 8,201 genes after removing unannotated transcripts and/or selecting one representative for highly correlated transcripts (correlation coefficient >0.7), then performed 10-fold cross-validation on six independent prediction methods using the entire cohort (Fig. 7 A and Table S7). The resulting predictive accuracy on cross-validation was first compared with models using only laboratory/clinical variables, then to models using laboratory/clinical variables combined with transcriptional data. Accuracy of transcriptional prediction ranged from 71.4% to 78.7% (mean accuracy: 74.2%), exceeding that of laboratory/clinical parameters (mean: 67.8%). When transcripts were combined with laboratory/clinical data, their mean accuracy was slightly higher (mean: 75.7%), suggesting that clinical parameters provided additive predictive value to transcriptional data. Receiver operating curves showed a similar trend (Fig. 7 B). The top 58 genes selected at least once in five of six models during cross-validation are depicted in Fig. 7 C. These include IFN-inducible and plasma cell–related transcripts, as well as transcripts related to platelet aggregation, angiogenesis, inflammation, and cell adhesion. Predicting OC using a similar machine learning approach yielded low accuracy, most likely due to the fact that OC includes a number of clinically heterogeneous complications and suggests that a larger sample size will be needed to develop independent classifiers for each of them. Overall, and although further validation in independent cohorts is warranted, our data suggest that early transcriptional changes in maternal blood might help predict PE in the context of SLE.

SLE is characterized by failed tolerance to self, and maintenance of a successful pregnancy requires unique tolerance of paternally encoded antigens. It is therefore not surprising that SLE pregnancy carries a risk of complications. Early biomarkers to predict pregnancy outcome and help inform treatment to decrease morbidity and mortality are needed. This study aimed to assess maternal blood transcriptome dynamics during healthy pregnancy, from embryo implantation to postpartum, and to determine molecular changes associated with outcome in the context of SLE pregnancies.

We found that healthy pregnancy is characterized by a significant modulation of major transcriptional networks associated with SLE pathogenesis, including IFN, plasma cell, lymphoid, myeloid, and erythroid cell–related pathways (Banchereau et al., 2016). Furthermore, these changes take place within days following embryo transfer, as detected in women undergoing ART. Importantly, immediately upon successful implantation an IFN-inducible transcriptionally spike is followed by down-regulation of the signature that persists through late pregnancy and PP. While this might be important to prevent the antiangiogenic milieu induced by IFNα (Andrade et al., 2015), down-regulation of the IFN response in human pregnancy could underlie the higher risk for certain infections such as influenza (Vanders et al., 2013), as also described in an allogeneic mouse pregnancy model (Engels et al., 2017).

Persistent down-regulation of the IFN signature compared with the increased SLE nonpregnant baseline is mirrored in noncomplicated SLE pregnancy and supports the physiological relevance of this finding, especially since the IFN signature remained more prominent in complicated lupus pregnancies. It also supports recent reports that activation of the type I IFN/IFNAR axis in pregnancy primes for inflammation-driven preterm birth in both mice and humans (Cappelletti et al., 2017; Engels et al., 2017). Whether the known association between SLE disease activity and pregnancy complications stems directly from increased IFN activity remains to be addressed (Molad et al., 2005; Clowse et al., 2008; Liu et al., 2012).

A recent study using mass cytometry revealed an increased STAT1 signaling response to ex vivo type I IFN in various innate and adaptive immune cells during pregnancy (Aghaeepour et al., 2017). Our results suggest, in turn, that increased ex vivo IFN-inducible signaling might reflect a lower threshold for activation due to in vivo transcriptional repression of this response during pregnancy. Indeed, monocytes from both human and mouse neonates have been reported to respond to TLR4 stimulation with strong induction of TIR-domain-containing adapter-inducing interferon-ß (TRIF)-dependent regulatory genes but weak activation of MyD88-dependent proinflammatory genes. This reciprocal shift was explained by the inverse baseline state of preactivated MyD88-dependent and barely detectable TRIF-dependent gene expression in neonatal versus adult monocytes (Ulas et al., 2017), which fits the transcriptional patterns we now describe in maternal blood during pregnancy.

In addition to the IFN signature, we detected sustained down-regulation of plasmablast/plasma cell signatures in both healthy and noncomplicated SLE pregnancies. These signatures, which are normally induced after infection and vaccination (Kardava et al., 2014; Kwissa et al., 2014; Athale et al., 2017), represent a robust biomarker of disease activity in pediatric SLE patients (Banchereau et al., 2016). Prolonged down-regulation of the plasma cell signature might be the result of either recruitment of antibody-secreting cells to the maternal–fetal interphase or active repression of terminal B cell differentiation programs. In fact, absolute counts of mature B cells are significantly lower, especially during the latest pregnancy stages (Kraus et al., 2012; Lima et al., 2016). This could be dependent on the inhibition of the IFN pathway, as type I IFN induces mature B cell survival and differentiation in multiple ways (Jego et al., 2003). Importantly, the plasma cell signature increased during the third trimester in SLE pregnancies associated with fetal complications. Thus, prospective studies to determine whether this signature represents an early biomarker of fetal complications are warranted.

Our study further reveals an up-regulation, in both healthy and noncomplicated SLE pregnancies, of innate immune signatures related to neutrophil, myeloid inflammation, and erythropoiesis pathways. Moreover, SLE pregnancy complicated by PE displayed the earliest (<15 wk) and most persistent up-regulation of the neutrophil signature, which correlated with an expansion of immature neutrophils. While glucocorticoid therapy has been linked to demargination of neutrophils, and this therapy was more frequently administered to SLE-PE patients at P1 compared with other groups (Table S6), daily glucocorticoid dose was low (prednisone <10 mg/d) and comparable between groups, and therefore unlikely to explain these changes (Guiducci et al., 2010).

The underlying mechanism and significance of neutrophil-related signatures in health and disease are not fully understood. These signatures might originate from more than one cell type, including myeloid-derived suppressor cells, which characteristically express high levels of Arginase-1 and inducible nitric oxide synthase (Monu and Frey, 2012). An expansion of immature neutrophils with myeloid-derived suppressor cell features has been reported in healthy pregnancy (Blazkova et al., 2017). Immature neutrophils are also expanded in nonpregnant SLE patients (Bennett et al., 2003; Blazkova et al., 2017). Furthermore, NETs associate with higher rates of endothelial dysfunction and thrombosis (Elkon and Wiedeman, 2012). NETs are also highly enriched in infiltrating placental intervillous spaces in non-SLE pregnancy complicated with PE and in SLE pregnancy compared with healthy controls (Marder et al., 2016). An in-depth elucidation of the complexity of the neutrophil lineage expansion during pregnancy, and of its value as a predictive marker of SLE-associated PE, will require further studies.

Increased transcription of genes linked to the erythroid lineage has not been reported before in pregnancy but has been detected in subsets of children with bacterial infections (Ardura et al., 2009), autoinflammatory diseases (Allantaz et al., 2007; Fall et al., 2007), and active SLE (Banchereau et al., 2016). A subset of CD71⁺ erythroid cells, which are enriched in neonatal mice and human cord blood, displays unique immunosuppressive properties through expression of Arginase-2 (Elahi et al., 2013). More recently, two distinct populations of tumor-inducible, immunosuppressive erythroblast-like cells were described in the context of advanced tumors. These cells facilitate tumor progression by secreting the neurotrophic factor artemin into the blood or through the production of reactive oxygen species (Han et al., 2018; Zhao et al., 2018). Whether similar cells and/or mechanisms play a role during pregnancy also needs to be further addressed.

Supporting previous reports, we observed increased T and NK cell signatures in P3 and P4 associated with both OC and PE. Indeed, an expansion of ICOS⁺ total CD4⁺ and Tfh cells as well as transitional B cells was detected at earlier time points by FACS. Tfh cell dysregulation contributes to autoimmunity and allograft rejection (Nankivell and Alexander, 2010; Liu et al., 2013; Ueno et al., 2015). Similar cells accumulate within the uterus and placenta in a mouse model of allogeneic pregnancy and expand significantly upon PD-L1 blockade, which leads to fetal resorption (Zeng et al., 2016).

Our study identified early transcriptional changes performing better than previously reported laboratory/clinical parameters to predict PE. Among them, PKD1 and Oculocerebrorenal syndrome of Lowe protein are involved in apoptosis and renal tubulogenesis (Boletta et al., 2000; Bastos et al., 2009), TLR5 activation inhibits attachment of human trophoblast cells to endometrial cells (Aboussahoud et al., 2010), NLRP12 regulates human preimplantation development (Tuncer et al., 2014), and LEFTY1 inhibits decidualization (Li et al., 2014) and TGFβ signaling (Ulloa and Tabibzadeh, 2001) and succinate receptor 1 (SUCNR1), which is linked to platelet activation (Macaulay et al., 2007) and to renovascular hypertension through activation of juxtaglomerular cells (Högberg et al., 2011; Mills and O’Neill, 2014). While this predictor list warrants validation, it uncovers potentially relevant mechanistic connections to SLE pregnancy complications. Recently, cell-free RNA transcripts in maternal blood have been used for predicting gestational age and preterm delivery (Ngo et al., 2018). These transcripts did not overlap with our cell-associated predictive transcripts, but the synergy of these noninvasive approaches should be explored in future prospective studies.

Our study identifies significant modulation of transcriptional and cellular immune pathways during healthy pregnancy that may facilitate implantation of the semi-allogeneic fetus. These pathways are also linked to SLE pathogenesis, and failure to modulate them is associated with lupus pregnancy complications. Our findings provide a framework for future studies aimed at investigating these immunological pathways and developing therapeutic strategies to improve health outcomes for mothers with SLE and their offspring.

PROMISSE is a multicenter, prospective observational study of pregnancies in women with SLE (four or more revised American College of Rheumatology criteria; Hochberg, 1997) and/or antiphospholipid antibodies (Miyakis et al., 2006), as well as healthy controls at low risk of adverse outcomes (one or more successful pregnancy; no prior fetal death; and fewer than two miscarriages at <10 WG; Buyon et al., 2015). All patients were recruited at <12 WG. This study focuses on patients enrolled between September 2003 and August 2013 from six sites in the United States and one in Canada. All mechanistic studies, including microarray and flow cytometry, were performed at the Baylor Institute for Immunology Research. The following institutional review committees from the recruiting organizations approved sample collection, sample receipt, and research performed as part of this study: Hospital for Special Surgery Institutional Review Board (IRB), Intermountain Health Care Urban Central Region IRB (Utah), NYU School of Medicine IRB, Oklahoma Medical Research Foundation IRB, The Johns Hopkins Medical Institutions (Western IRB), The University of Chicago IRB, Mt. Sinai Hospital’s Research Ethics Committee (Toronto, Canada), and University of Utah IRB. Patients with SLE had inactive or mild/moderate activity at screening (mean SLEPDAI score = 3.01). Exclusion criteria were prednisone >20 mg/d, urinary protein/creatinine ratio >1,000 mg/g, erythrocyte casts, serum creatinine level >1.2 mg/deciliter, diabetes mellitus, and blood pressure >140/90 mm Hg at screening. Detailed medical and obstetrical information and serial blood samples were collected each month during the course of pregnancy. IRBs approved the protocol, and written, informed consent was obtained from all participants.

Adverse pregnancy outcomes were defined as one or more of the following: (1) PE at any time (ACOG Committee on Practice Bulletins–Obstetrics, 2002); (2) fetal death at >12 WG unexplained by chromosomal abnormality, anatomical malformation, or congenital infection; (3) neonatal death before hospital discharge due to complications of prematurity; (4) preterm delivery or termination of pregnancy at <36 WG due to growth restriction or placental insufficiency; and (5) small for gestational age of less than the fifth percentile at birth. We subdivided patients into those with PE and those with OC.

Our target sample size was 40 SLE pregnancies without complications, 40 SLE pregnancies with complications (20 with PE; 20 with OC), and 40 healthy control subjects who had complete sets of PreAnalytiX blood RNA (PAXgene) tubes (Qiagen) for each time point of interest. However, to avoid dropouts due to sample processing, a few extra patients were included in each group. We identified 46 SLE patients with complications who had complete sets of PAX gene tubes available for analysis and collected at the following intervals: <16 WG, 16–23 WG, 24–31 WG, 32–40 WG (if still pregnant), and 8–20 wk PP (Fig. 1 B). For comparison groups, we selected the first 46 SLE-NC and 43 H-P. We also enrolled 23 H-NP controls. Laboratory values and demographics of SLE patients are summarized in Table S6.

For SLE-NP and additional H-NP, 82 subjects (62 SLE-NP and 20 H-NP individuals) were obtained from NCBI GEO (accession no. GSE49454). These samples were originally hybridized on the same platform (Illumina HT-12 V4 beadchips) and in the same facility. Exclusion criteria included (1) current treatment with mycophenolate mofetil or cyclophosphamide, (2) active renal flare, (3) patient >43 yr, and (4) male patient. When patients had multiple visits, the sample with the lowest SLEPDAI score was chosen, resulting in the selection of 20 unique SLE-NP and 11 H-NP samples. There were no demographic, clinical, or disease activity differences between pregnant and nonpregnant patient groups (Table S6).

An independent cohort of healthy women undergoing ART was recruited at the University of Connecticut Center for Advanced Reproductive Services from November 2014 to July 2015. The protocol was approved by the Independent Review Board of the University of Connecticut Health Center (IRB no.15-057J-3). Subjects who were <37 yr of age with normal ovarian reserve and undergoing blastocyst transfer were included in the study. Women with recent vaccination (within 14 d), chronic infection, or autoimmune disorders were excluded. Whole blood was collected in Tempus tubes on the day of embryo transfer, midluteal, and on the day of scheduled pregnancy test. Those who conceived (defined as positive pregnancy test) were asked to submit an additional sample between 6 and 8 WG (P0). Of the 40 subjects enrolled, 25 were included in the final analysis. The clinical protocol from seven subjects was terminated before embryo transfer, and samples from eight subjects failed quality control.

Whole-blood samples were collected throughout pregnancy from 50 individuals from PROMISSE for flow cytometry. This convenience sample, based on practical conditions of patient availability, yielded 125 blood samples: P1–4 and PP for H-P (n = 11) and SLE-NC (n = 30) and, due to scarcity and early delivery in complicated group, P1–2 for SLE-PE (n = 5) and SLE-OC (n = 4) samples were collected. Acid citrate dextrose–anticoagulated whole blood was aliquoted at 200 μl/tube to six 5-ml Falcon tubes (BD Biosciences) for control staining and phenotyping of B cells, monocytes and granulocytes, Tfh cells, dendritic cells, T regulatory cells, and NK cells. Blood was incubated with monoclonal antibodies for 15 min at room temperature, followed by RBC lysis with BD FACSLyse Buffer (BD Biosciences) for 10 min at room temperature in the dark, and two washing steps with PBS. Samples were acquired on a custom BD LSR Fortessa. The gating strategy used herein was previously described (Maecker et al., 2012; Obermoser et al., 2013). Data were analyzed in Flowjo 9.8 software. FACS antibodies are listed in Table S8.

Total RNA was extracted as previously described (Berry et al., 2010). Amplified RNA was hybridized to Illumina HT-12 V4 beadchips (47,231 probes) and scanned on an Illumina Beadstation 500. Signal-intensity values were generated using Illumina’s GenomeStudio version 2011.1 with the Gene Expression Module v1.9.0 (Illumina). Background subtraction and average normalization were performed using GenomeStudio. The microarray dataset described in this manuscript is deposited in the NCBI GEO (accession no. GSE108497).

RNA was isolated and globin-reduced using the Tempus Spin RNA kit (Applied Biosystems/Ambion) and the GLOBINclear kit (Applied Biosystems/Ambion), respectively. Illumina barcoded stranded RNA-seq libraries were prepared and sequenced using either the Illumina HiSeq2500 or NextSeq500 targeting 30 million paired-end reads (2 × 75 bp) per sample. FASTQ files were created using the Illumina CASAVA v1.8.2 pipeline. The sequencing reads generated by the pipeline were further subjected to quality control using NGSQC toolkit v2.3 and samples with at least 70% of the nucleotides having base quality >30 were included in the analysis. High-quality reads were mapped to human transcriptome from ENSEMBL build 70 based on the NCBI human genome version 37. Alignments were made with TopHat v2.0 and then subsequent transcript and gene-level counts, fragments per kilobase of transcript per million reads, and transcripts per million values were obtained using RSEM using the default parameters. The RNA-seq data for ART is deposited in BioProject (accession no. PRJNA427177).

Two linear mixed models were developed to identify DETs in microarray (A) or differentially expressed flow subsets (B). To exploit longitudinal and repeated measures, all models consisted of a patient-specific random intercept and specific fixed effects. Spatial power covariance structures were used for all linear mixed models to explain regular intervals between repeated measurements. Estimate statements were used to test for comparisons of interest. Models used the SAS MIXED procedure using sp(pow) type and were run in JMP Genomics 6.0 (SAS Institute, Inc.). Detailed SAS codes are provided in Table S10.

The first model was designed to identify DETs correlated with disease, pregnancy, complication group, and time point. APO = Adverse Pregnancy Outcome (PE + OC).

The second linear mixed model was developed to identify FACS subsets correlated with complication group and tim point.

Our dataset was normalized with batch correction method using JMP Genomics 6.0 after standardizing for batch (SAS Institute, Inc.). Briefly, datasets from two batches were normalized by establishing a batch profile based on averaging across within-batch-level H-NP controls and then applying this profile to correct values across batches to all the other samples. K-mean clustering is used for grouping batch profiles into clusters.

A preexisting framework of 260 coexpressed transcriptional modules was used (Chaussabel and Baldwin, 2014). For simplicity, only the 97 modules forming the first seven rounds of selection of the module set were represented. For each module, the percentage of probes significantly up- or down-regulated in the pregnant group was calculated based on the linear mixed-model result. The module score was calculated based on the percentage of genes within a module that are significantly (FDR < 0.05) upregulated minus the percentage of genes within a module that are significantly down-regulated (within a given comparison). Gene lists for Baylor whole blood modules are listed in Table S9.

Q-Gen is a generalization of the QuSAGE algorithm (Yaari et al., 2013) that integrates statistics from mixed-model analysis by conducting necessary estimation adjustments (Turner et al., 2015). Modules with at least 15 transcripts were considered for QuSAGE analysis. Significantly enriched gene sets (blood modules) were selected based on two criteria: (1) FDR < 0.05 and (2) FC ≥ 1.3 or ≤ 0.7 in the comparison considered in the training set.

Predictive modeling was performed using the Caret package. ROC curves are generated using the pROC.

GenBart files are generated using genBART R package. BART (Biostatistical Analysis Reporting Tool) is a user-friendly, point-and-click, R shiny application. The analysis results were formatted and combined to an R datafile to be uploaded and efficiently examined (available at https://data.mendeley.com/datasets/hg596dwdr6/1; Hong, 2018). BART provides users the ability to easily view, modify, and download the tables and figures generated by the app.

Fig. S1 shows line graphs of significantly differentiated modules in SLE-complicated groups in Fig. 6 B normalized to H-NP. Fig. S2 compares blood module signatures of SLE-PE and SLE-OC with or without adjustment by SLEPDAI score at initial visit. Fig. S3 represents a summary of the main features of genBart shiny app. Table S1 shows demographics, embryo transfer method, and estradiol level before embryo transfer categorized by pregnancy status for the ART cohort. Table S2 shows Q-Gen results based on a mixed-model analysis using blood modules as gene sets. Table S3 shows absolute counts of 84 blood cell subsets measured by flow cytometry. Table S4 shows flow cytometry populations (n = 43) significantly changed in a linear mixed model (P ≤ 0.05). Table S5 shows Q-Gen results based on a linear mixed model using blood modules as gene sets for the ART cohort. Table S6 shows demographics, laboratory values, SLE activity, and medications of SLE patients from the PROMISSE cohort during pregnancy. Table S7 lists predictors’ variable importance and their significant index of 8,201 genes in five PE prediction models. Table S8 shows an antibodies list used for flow cytometry. Table S9 shows gene lists corresponding to whole-blood modules. Table S10 shows detailed SAS codes.

We are grateful to the SLE patients and healthy subjects who participated in our study and M.D. Lockshin for help with recruiting patients.

This work was supported by the National Institute of Arthritis and Musculoskeletal and Skin Diseases/National Institutes of Health awards R01 AR49772 and R01 AR49772-07S2 (J.E. Salmon), R01 AR69572 and AR43727 (M. Petri), National Institute of Allergy and Infectious Diseases/National Institutes of Health award U19 AI082715 (V. Pascual), P50AR070594 (V. Pascual and J. Banchereau), the Morris and Alma Schapiro Fund (J.E. Salmon), the Baylor-Scott & White Health Care Research Foundation, and the Kathryn W. Davis gift to the Jackson Laboratory.

Author contributions: S. Hong and R. Banchereau analyzed data and wrote the manuscript. B-S.L. Maslow, J. Nulsen, and J. Banchereau recruited and generated the ART data. S. Hong, B-S.L. Maslow, and D. Nehar-Belaid analyzed the ART data. M.M. Guerra collected clinical data and revised the manuscript. J. Baisch, D.W. Branch, T.F. Porter, A. Sawitzke, C.A. Laskin, J.P. Buyon, J.T. Merrill, L.R. Sammaritano, and M. Petri followed PROMISSE patients and collected clinical data and samples. S. Hong, J. Turner, J. Cardenas, and D. Blankenship conducted statistical analysis, developed models, and developed the GenBart interface. E. Gatewood, A-M. Cepika, M. Ohouo, T.W. Kim, and G. Obermoser generated and interpreted FACS analysis. E. Anguiano supervised microarray. V. Pascual, J.E. Salmon, and J. Banchereau supervised the study, interpreted data, and wrote the manuscript.

The authors declare no competing financial interests.