Gene expression responses to Zika virus infection in peripheral blood mononuclear cells from pregnant and non-pregnant women.

Congenital Zika syndrome is caused by mother-to-fetus transmission of the Zika virus (ZIKV). Peripheral blood mononuclear cells (PBMCs) are permissive to ZIKV infection and may carry ZIKV to the placenta. To identify pregnancy-related differences in PBMC responses against ZIKV infection, we compared gene expression profiles of ZIKV-infected and non-infected PBMCs cultured from pregnant and non-pregnant women. ZIKV-infected pregnant conditions generally overexpressed M1-shifted pro-inflammatory responses and underexpressed M2-shifted anti-inflammatory responses. Additionally, transcripts involved in osteoclast differentiation and cardiac myopathies were upregulated following ZIKV infection. Our results suggest potential roles of pregnancy-induced immune dysregulation in shaping neonatal pathology associated with ZIKV infection.

INTRODUCTION

Pregnancy is a unique immunological state that generally protects the fetus, but susceptibility during pregnancy to various infections has been observed. In parasitic infections, such as malaria and toxoplasmosis, increased host susceptibility and microbial virulence during pregnancy have been reported (Piao et al., 2018; Tuikue Ndam et al., 2018). Pregnancy is also associated with increased susceptibility to viral infections, including influenza virus, hepatitis E virus, herpes simplex virus, and Zika virus (ZIKV) (Kourtis et al., 2014; Silasi et al., 2015). Infectious mechanisms vary, and some involve the placenta and membranes while others could impact the fetus through immunological perturbations.

ZIKV is an epidemic‐causing pathogen transmissible through Aedes mosquitoes, sexual contact, and transplacental transfer (Paz‐Bailey et al., 2019). Transplacental ZIKV transmission causes congenital Zika syndrome characterized by birth anomalies, including microcephaly, brain damage, ocular defects, and neuromuscular impairments (Mawson, 2016). ZIKV targets peripheral blood mononuclear cells (PBMCs), particularly monocytes and, to a lesser extent, dendritic cells (Foo et al., 2017; Michlmayr et al., 2017). These circulating cells may facilitate ZIKV delivery to the placenta (Michlmayr et al., 2017) and affect placental development through the regulation of pro‐inflammatory M1 macrophages and anti‐inflammatory M2 macrophages ratio (Foo et al., 2017). PBMCs and monocytes from pregnant women produced aggravated M1‐shifted pro‐inflammatory responses when experimentally infected with the African‐lineage ZIKV strain MR766, but M2‐shifted immunosuppression responses when infected with the Asian‐lineage ZIKV strain H/PF/2013 (Foo et al., 2017). This suggests that pregnant women are more susceptible to the ZIKV strain H/PF/2013, although the generalizability of this finding across other ZIKV strains has not yet been explored.

In this study, we investigated whether PBMCs from pregnant women show differential immune responses to ZIKV strain Puerto Rico PRVABC59 of Asian lineage compared to PBMCs from non‐pregnant women. Using PBMCs isolated from uninfected pregnant or non‐pregnant women, we cultivated and experimentally infected them with the ZIKV strain. Using RNAseq, we evaluated gene expression changes pre‐ and post‐ZIKV infection in vitro using PBMCs derived from pregnant and non‐pregnant women, followed by a comparison of gene expression profiles between ZIKV‐infected pregnant conditions and ZIKV‐infected non‐pregnant conditions. With the sequenced data, we also explored whether other non‐immune pathways potentially relevant to pregnancy and birth outcomes were differentially expressed across ZIKV infection, pregnancy, and/or IL‐3 stimulation conditions.

MATERIALS AND METHODS

PBMC isolation and infection with ZIKV virus

PBMCs were isolated by Ficoll‐Paque density gradient purification using 10 ml of fresh blood from two healthy donors who were seven and 15 weeks pregnant and two non‐pregnant healthy female donors who had no underlying conditions and who provided written informed consent (Table A1). Research procedures adhered to the Declaration of Helsinki and were approved by the university Institutional Review Board (USF IRB protocol Pro00029803). PBMCs were seeded in 24‐well plates at a density of 1 million cells per well in 1 ml of DMEM containing 5% fetal bovine serum (FBS) with penicillin/streptomycin and amphotericin B. PBMCs from each donor were sub‐cultured in triplicates under the following conditions: PBMC control, PBMC + ZIKV, PBMC + interleukin (IL)‐3 (30 µM), and PBMC + ZIKV + IL‐3. Interleukin‐3, a hematopoietic growth factor, was added to stimulate the survival, proliferation, and differentiation of hematopoietic cell types from myeloid progenitor cells (Stocking & Ostertag, 1990) throughout the cell culture period (48 h). After 5 h of incubation at 37°C, PBMCs were infected with 150 µl of ZIKV at a multiplicity of infection (MOI) of 0.5 (with MOI = pfu of ZIKV per number of PBMCs), followed by incubation at 37°C for 48 h. Light microscopy detected no cell lysis within 48 h.

Table A1

Participants' information

Age	Ethnicity	Race	Birth country	Weeks pregnancy	Pregnancy status
25	Non‐Hispanic	White	USA	NA	Non‐pregnant
21	Hispanic	White	Colombia	NA	Non‐pregnant
38	Non‐Hispanic	Black	USA	7	Pregnant
38	Hispanic	White	Venezuela	15	Pregnant

ZIKV virus culture and purification

ZIKV Puerto Rico strain PRVABC59 (ATCC VR‐1843) was propagated in Vero E6 cells (ATCC CRL‐1586) grown in DMEM medium supplemented with 5% FBS. Vero E6 cells at 80% confluency in a T25 flask were washed for 10 min with 1 ml of serum‐free DMEM before infection with 2 ml of ZIKV in DMEM at MOI = 0.5 for 1 h at 37°C. Subsequently, 5 ml of DMEM with 5% FBS was added and cultures were incubated at 37°C and 5% CO2 atmosphere. ZIKV was collected 5–10 days post‐infection when ~80% of host cells were lysed. The supernatant containing ZIKV particles was aspirated, mixed with 700 µl of 10× sucrose phosphate glutamate (SPG) virus freezing medium (10× SPG: 2.18 M sucrose, 38 mM KH2PO4/72 mM K2HPO4, 49 mM l‐glutamate), and stored in aliquots at −80°C. ZIKV infectious titer was determined in 24‐well plates using the Tissue Culture Infective Dose resulting in 50% of lysed Vero E6 cells (TCID50) endpoint dilution assay protocol at 5 days endpoint post‐infection and pfu were calculated by the formula pfu/ml = 0.5 × TCID50 (Hierholzer & Killington, 1996).

RNA extraction, sequencing, and data analyses

Forty‐eight hours post‐infection, the PBMC cell lysate from each well was homogenized in Qiagen's RLT buffer with 1% beta‐mercaptoethanol and centrifuged in the QIAshredder spin column for 2 min. On‐column DNA digestion was performed using Qiagen's RNase‐free DNase I before total RNA extraction using Qiagen's RNeasy Mini kit. Extracted RNA showed 2.07–2.14 OD260/OD280 ratio and 8.9–9.5 RNA integrity number. For each culturing condition, 1.6–2.5 µg of total RNA per sample were pooled and sent to Novogene (Chula Vista, CA). RNA quality was assessed using NanoDrop spectrophotometer, agarose gel electrophoresis, and Agilent 2100 bioanalyzer before library preparation using the NEBNext® UltraTM RNA library preparation kit and sequencing on Illumina's NovaSeq 6000 2 × 150 bp platform. As such, each sequenced library comprises RNA pooled from all pregnant or non‐pregnant subjects from each culturing condition.

Novogene's analysis pipeline removed reads containing adaptors, >10% undetermined bases, and ≤5 quality score of >50% bases (Table A2). The remaining reads were mapped to the Ensembl (Cunningham et al., 2019) reference human assembly GRCh37 (hg19) using Spliced Transcripts Alignment to a Reference (STAR) v2.5 (Dobin et al., 2013) (Table A3). Sequencing depth was calculated using SAMtools v1.4 (Li et al., 2009) (Table A2). The number of reads mapped to each gene was counted using HTSeq v0.6.1 (Anders et al., 2015) and normalized to transcript per million (TPM) to identify the most abundantly expressed transcripts across all libraries. TPM values were also used to deconvolute PBMC transcriptomic data computationally into specific immune cell types using ABsolute Immune Signal (ABIS) (Monaco et al., 2019). Principal component analysis (PCA) was performed on the ABIS output using R v3.6.1.

Table A2

Summary of read qualities before and after trimming

Library	# Raw reads	# Clean reads	# Bases (raw reads)	# Bases (clean reads)	Q20 (%)	Q30 (%)	Average depth	GC content (%)
NP_Ctr	74,395,980	70,713,678	11.2G	10.6G	97.2	92.7	7.1	49.2
NP_ZV	65,876,434	63,735,562	9.9G	9.6G	97.2	92.7	6.3	50.3
NP_IL3	62,863,948	56,729,120	9.4G	8.5G	97.1	92.6	5.5	50.6
NP_ILZV	50,445,536	49,018,202	7.6G	7.4G	97.3	92.9	4.9	50.4
P_Ctr	57,073,364	55,022,010	8.6G	8.3G	97.3	93.1	5.4	49.9
P_ZV	54,009,684	52,060,670	8.1G	7.8G	97.0	92.3	5.1	50.5
P_IL3	53,231,182	51,597,452	8.0G	7.7G	97.2	92.7	5.2	49.8
P_ILZV	57,493,344	55,393,380	8.6G	8.3G	97.1	92.6	5.5	50.5

Table A3

Summary of read‐mapping results

Library	Total reads	Total mapped	Multiple mapped	Uniquely mapped	Reads mapped to +strand	Reads mapped to −strand	Reads mapped to splice junctions
NP_Ctr	70,713,678	65,256,562 (92.3%)	1,408,190 (2.0%)	63,848,372 (90.3%)	31,924,186 (45.2%)	31,924,186 (45.2%)	21,119,672 (29.9%)
NP_ZV	63,735,562	58,163,546 (91.3%)	1,155,376 (1.8%)	57,008,170 (89.4%)	28,504,085 (44.7%)	28,504,085 (44.7%)	21,073,355 (33.1%)
NP_IL3	56,729,120	51,078,268 (90.0%)	1,111,686 (2.0%)	49,966,582 (88.1%)	24,983,291 (44.0%)	24,983,291 (44.0%)	16,753,944 (29.5%)
NP_ILZV	49,018,202	44,714,808 (91.2%)	887,188 (1.8%)	43,827,620 (89.4%)	21,913,810 (44.7%)	21,913,810 (44.7%)	15,728,909 (32.1%)
P_Ctr	55,022,010	48,832,694 (88.8%)	1,111,032 (2.0%)	47,721,662 (86.7%)	23,860,831 (43.4%)	23,860,831 (43.4%)	13,766,361 (25.0%)
P_ZV	52,060,670	46,541,166 (89.4%)	937,810 (1.8%)	45,603,356 (87.6%)	22,801,678 (43.8%)	22,801,678 (43.8%)	16,797,986 (32.3%)
P_IL3	51,597,452	47,114,396 (91.3%)	1,037,214 (2.0%)	46,077,182 (89.3%)	23,038,591 (44.7%)	23,038,591 (44.7%)	15,564,627 (30.2%)
P_ILZV	55,393,380	49,852,760 (90%)	1,007,890 (1.8%)	48,844,870 (88.2%)	24,422,435 (44.1%)	24,422,435 (44.1%)	18,672,101 (33.7%)

For differential expression analyses, raw counts were imported to edgeR v3.26.5 (Robinson et al., 2010) and filtered to retain 18,594 genes with >0.5 count per million in at least two libraries. The filtered count data across libraries were normalized by the trimmed mean of M values method to adjust for variations in library sizes (depths) and composition (Robinson & Oshlack, 2010). Pearson correlation analysis was performed on fragments per kilobase million normalized counts of each library. Multi‐dimensional scaling analysis was performed by converting the gene count matrix of each library to log counts per million. For each pairwise comparison between libraries, the top 500 genes that distinguish the libraries were used to calculate the Euclidean distances (root‐mean‐square deviation) between each pair of samples. Differentially expressed genes were predicted using the likelihood ratio test with thresholds of p ≤ 0.05 adjusted for false discovery rate (FDR), based on fitted negative binomial generalized linear models. The common dispersion of the main model was estimated using edgeR's estimateGLMCommDisp function with the parameters: method = “deviance,” robust = TRUE, and subset = NULL. This is recommended in the edgeR user's guide for our reduced design matrix, where libraries P_ZV and P_ILZV were treated as replicates and libraries NP_ZV and NP_ILZV were treated as replicates in the absence of biological replication. Statistical power was calculated using the RNASeqPower R package (https://bioconductor.org/packages/release/bioc/html/RNASeqPower.html). Differentially expressed genes and pathways were annotated using UniProt (UniProt, 2019), GeneCard (Stelzer et al., 2016), and the KEGG mapper tool (Kanehisa & Sato, 2019). Classical M1 macrophage and M2 macrophage marker genes were compiled from the literature (Chavez‐Galan et al., 2015; Foo et al., 2017; Roszer, 2015; Sica & Mantovani, 2012).

RESULTS

Predicted immune cell proportions from PBMC transcriptomic data

We analyzed eight PBMC transcriptomes from total RNA pooled from pregnant (P) and non‐pregnant (NP) women‐derived PBMCs without ZIKV or IL‐3 (P_Ctr and NP_Ctr), with ZIKV only (P_ZV and NP_ ZV), with IL‐3 only (P_IL3 and NP_IL3), and with IL‐3 and ZIKV (P_ILZV and NP_ILZV). Each transcriptomic library consisted of ~57 ± 7 million reads (~6 ± 1× sequencing depth) on average after quality filtering (Table A1). An average of 91 ± 0.01% of these clean reads mapped to the reference human genome, and ≤2% of the reads mapped to more than one location in the genome (Table A2).

From the transcriptomic data, we estimated immune cell proportions using the ABsolute Immune Signal (ABIS) deconvolution algorithm (Monaco et al., 2019). Higher β coefficients, corresponding to higher proportions and mRNA abundances (Monaco et al., 2019), were observed for CD4 memory T cells in uninfected libraries (control and IL3‐only) relative to ZIKV+ samples with and without IL‐3 (Figure 1a). In both uninfected and infected conditions, libraries from pregnant women showed higher β coefficients for CD8 naïve T cells and CD4 memory T cells, but lower β coefficients for CD4 naïve T cells, compared to those from non‐pregnant women (Figure 1a). In pregnant women, β coefficients for CD8 naïve T cells and CD4 memory T cells were higher, but β coefficients for CD4 naïve T cells were lower in uninfected conditions relative to ZIKV‐infected conditions (Figure 1a). PCA analysis of the β coefficients identified CD8 naïve T cells, CD4 memory T cells, and CD4 naïve T cells as the strongest contributors to the observed variations between pregnant and non‐pregnant conditions (Figure 1b). The PCA analysis also grouped libraries P_ZV and P_ILZV into one cluster and libraries NP_ZV and NP_ILZV into another cluster, suggesting relatively smaller effects of IL‐3 stimulation on predicted immune cell types under ZIKV‐infected conditions (Figure 1b). Contrarily, no clear clustering pattern was observed for immune cell proportions of non‐infected cultures (Figure 1b). Immune cell proportions that varied the most in ZIKV‐infected PBMCs with and without IL‐3 treatment included CD8 naïve T cells (8.1% decrease in NP_ILZV compared to NP_ZV and 17.0% decrease in P_ILZV compared to P_ZV) and CD4 memory T cells (7.8% decrease in NP_ILZV compared to NP_ZV and 13.2% decrease in P_ILZV compared to P_ZV).

FIGURE 1

(a) Deconvoluted proportions of immune cell types from PBMC transcriptomic data. The β‐coefficient is an indicator of the proportion and mRNA abundance of a specific immune cell type (Monaco et al., 2019). LD, low‐density granulocytes; MAIT, mucosal‐associated invariant T cells; mDCs, myeloid dendritic cells; monocytes C, classical monocytes; monocytes NC+I, non‐classical and intermediate monocytes; pDCs, plasmacytoid dendritic cells; T gd non‐Vd2, gamma delta T cell without T‐cell receptor Vδ2 antigen; T gd Vd2, gamma delta T cell with T‐cell receptor Vδ2 antigen. (b) Principal component analysis (PCA) biplot showing the effects of predicted immune cell types on sample clustering, based on β‐coefficients that estimate cell proportions and mRNA abundances

Transcriptome profiles across libraries

The most abundant transcripts across all libraries participate in protein biosynthesis, respiratory, cytoskeletal, and immune functions (Figure A1). Multi‐dimensional scaling (MDS) analysis based on the root‐mean‐square log2‐fold change for the top 500 genes that distinguish each library pair showed similar clustering patterns to the PCA analysis on predicted immune cell proportions (Figure 2), although we note that both analyses were derived from the same gene count data. In addition to their clustering on the MDS plot (Figure 2), libraries P_ZV and P_ILZV and libraries NP_ZV and NP_ILZV shared >99% pairwise Pearson correlations between their transcript counts. Because of their gene expression profile similarities, for all downstream analyses, we treated libraries P_ZV and P_ILZV as replicates and libraries NP_ZV and NP_ILZV as replicates (Figure 2). Based on our experimental design (samples = 2, depth = 6, and estimated biological coefficient of variation of the fitted model (CV = 0.072)), our study was able to detect a 2‐fold change in gene expression between categories with 0.39 statistical power at α = 0.05.

FIGURE A1

Log10‐transformed transcript per million (TPM) of gene products of the 35 most abundantly expressed protein‐coding transcripts across all sequenced samples

FIGURE 2

(a) Heat map of pairwise Pearson correlations, based on fragments per kilobase million, across sequenced libraries; (b) multi‐dimensional scaling plot showing distances between libraries, based on the root‐mean‐square log2‐fold change for the top 500 genes that distinguish each library pair. Samples treated as replicates for differential expression analyses are boxed in (a) and circled in (b)

Differentially expressed immune‐related transcripts before and after ZIKV infection

We first sought to identify differentially expressed immune‐related transcripts following ZIKV infection in pregnant (P_ZV and P_ILZV compared to P_Ctr) and non‐pregnant conditions (NP_ZV and NP_ILZV compared to NP_Ctr), respectively. In pregnant conditions, 4747 genes were upregulated, and 3358 genes were downregulated with ZIKV infection. In non‐pregnant conditions, 3853 genes were upregulated, and 2140 genes were downregulated with ZIKV infection. Among these, 3278 genes were commonly upregulated, and 1322 genes were commonly downregulated with ZIKV infection in pregnant and non‐pregnant conditions. Commonly overexpressed immune‐related transcripts encoded for C‐C motif chemokines (CCL1‐CCL4, CCL7, CCL8, CCL18, CCL20, and CCL24), C‐C chemokine receptors (CCR2‐CCR3), C‐C chemokine receptor‐like 2 (CCRL2), C‐X‐C motif chemokines (CXCL1, CXCL2, CXCL5, CXCL9‐CXCL11, and CXCL16), interleukins (IL1A, IL1B, IL8, IL10, and IL18), interleukin‐1 receptor accessory protein (IL1RAP), toll‐like receptors (TLR2 and TLR4), and tumor necrosis factor receptor superfamily members (TNFRSF1A, TNFSF10, and TNFSF14). Interferon receptor genes, interferon‐inducible genes, interferon regulatory genes, and interferon signaling genes (Foo et al., 2017) were also commonly upregulated with ZIKV infection in both pregnant and non‐pregnant conditions. These transcripts encoded interferon gamma receptors 1 (INFGR1 and INFGR2), interferon‐alpha/beta receptor 1 (IFNAR1), interferon‐alpha‐inducible proteins (IFI16, IFI27, IFI30, IFI35, IFI44, and IFI6), interferon‐alpha‐inducible protein‐like proteins (IFI27L2 and IFI44L), interferon‐induced proteins with tetratricopeptide repeats (IFIT1‐IFIT3 and IFIT5), interferon‐induced transmembrane proteins (IFITM1‐IFITM3), interferon‐induced helicase C domain‐containing protein 1 (IFIH1), interferon‐induced GTP‐binding proteins (MX1 and MX2), interferon‐induced antiviral 2′‐5′‐oligoadenylate synthases (OAS1‐OAS3), 2′‐5′‐oligoadenylate synthase‐like protein (OAL), interferon‐induced antiviral protein (RSAD2), interferon‐related developmental regulator 1 (IFRD1), interferon regulatory factors (IRF1, IRF5‐IRF8), Janus kinase (JAK2), Janus kinase and microtubule‐interacting protein 2 (JAKMIP2), toll‐like receptor‐interacting protein (MYD88), and nuclear factor kappa B subunit p50 (NFKB1). Immune‐related transcripts commonly downregulated with ZIKV infection included CCR4 and IL16.

Pro‐inflammatory M1 macrophage‐related and anti‐inflammatory M2 macrophage‐related marker genes (Chavez‐Galan et al., 2015; Foo et al., 2017; Roszer, 2015; Sica & Mantovani, 2012) were also differentially expressed following ZIKV infection in pregnant and non‐pregnant conditions (Figure A2). Transcripts overexpressed in pregnant conditions, but not in non‐pregnant conditions, following ZIKV infection were M1‐associated (interferon gamma [IFNγ]; logFC = 1.9 and CCR7; logFC = 0.2), M1/M2‐associated (tumor necrosis factor [TNF]; logFC = 0.5 and CCL5; logFC = 0.3), and M2‐associated (CCL17; logFC = 3.0; Figure A2). No M1‐related transcript was overexpressed only in non‐pregnant conditions following ZIKV infection. Transcripts upregulated following ZIKV infection in non‐pregnant conditions, but not in pregnant conditions, were M1/M2‐associated (IL6; logFC = 2.3 and HLA‐DRB5 encoding DR beta 5 chain of HLA class II histocompatibility antigen; logFC = 1.1) and M2‐associated (VEGFA encoding vascular endothelial growth factor A; logFC = 0.9; Figure A2). In pregnant conditions, one M2‐related transcript encoding cell surface glycoprotein CD200 receptor 1 (CD200R1; logFC = −0.5) and another M1/M2‐related transcript encoding mitochondrial arginase‐2 (ARG2; logFC = −0.7) were downregulated in ZIKV+ PBMCs compared to non‐infected PBMCs (Figure A2). Transcripts downregulated in ZIKV+ PBMCs compared to non‐infected PBMCs in non‐pregnant condition were exclusively M2‐related and included arginase‐1 (ARG1; logFC = −2.2), CD200R1 (logFC = −0.8), C‐C motif chemokine 22 (CCL22; logFC = −0.5), and the interferon regulatory factor 4 (IRF4; logFC = −0.3; Figure A2). Many differentially expressed immune‐related genes overlap with those identified in a separate analysis comparing PBMC gene expression between infected (P_ZV, NP_ZV, P_ILZV, and NP_ILZV) and non‐infected (P_Ctr, NP_Ctr, P_IL3, and NP_IL3) cultures, where pregnant, non‐pregnant and IL‐3+/− libraries were combined (Appendix 1, Figure A3, and Figure A4). This implies that these immune‐related genes were indeed commonly differentially expressed following ZIKV infection and that their differential expression was likely not IL‐3‐dependent.

FIGURE A2

Log2‐fold changes (logFC) of differentially expressed transcripts associated with M1‐shifted pro‐inflammatory responses and M2‐shifted anti‐inflammatory responses between ZIKV+ and ZIKV− PBMCs from (a) pregnant and (b) non‐pregnant women. INHBA, inhibin beta A chain; MRC1, macrophage mannose receptor 1; CCL, C‐C motif chemokine; IDO1, indoleamine 2,3‐dioxygenase 1; CLEC10A, C‐type lectin domain family 10 member A; TGM2, protein‐glutamine gamma‐glutamyltransferase 2; IL1RN, interleukin‐1 receptor antagonist protein; EPAS1, endothelial PAS domain‐containing protein 1; CD, leukocyte differentiation antigen (cluster of differentiation); FCGR1A, high affinity immunoglobulin gamma Fc receptor I; NR1H3, oxysterols receptor LXR‐alpha; TLR, toll‐like receptor; IL, interleukin; ILR, interleukin receptor; CSF, granulocyte‐macrophage colony‐stimulating factor; P2RY2, P2Y purinoceptor 2; IFNG, interferon gamma; SOCS1, suppressor of cytokine signaling 1; STAT, signal transducer and activator of transcription; BTK, Bruton's tyrosine protein kinase; IRF, interferon regulatory factor; HLA, human leukocyte antigen; IL1R1, interleukin‐1 receptor type 1; HIF1A, hypoxia‐inducible factor 1‐alpha; PTGS2, prostaglandin G/H synthase 2; NFKB1, nuclear factor NF‐kappa‐B p105 subunit; CCR, C‐C chemokine receptor; TNF, tumor necrosis factor; TGFB1, transforming growth factor beta‐1 proprotein; CD200R1, cell surface glycoprotein CD200 receptor 1; ARG, arginase; VEGFA, vascular endothelial growth factor A

FIGURE A3

Log2‐fold changes (logFC) of the 35 most differentially expressed transcripts between ZIKV+ (P_ZV, NP_ZV, P_ILZV, and NP_ILZV) and ZIKV− (P_Ctr, NP_Ctr, P_IL3, and NP_IL3) PBMCs. Solid bars represent regular gene transcripts, cross‐hatched bars represent pseudogenes and stippled bars long non‐coding RNAs. Abbreviations for regular gene transcripts: TRAJ4, T‐cell receptor alpha joining; CYP453; cytochrome P450 4F3; GPR52, G‐protein‐coupled receptor 52; OGN, mimecan; TAS2R19, taste receptor type 2 member 19; SIGLEC1, sialoadhesin; PDGFRL, platelet‐derived growth factor receptor‐like protein; IFIT, interferon‐induced protein with tetratricopeptide repeats; APOBEC3A, apolipoprotein B mRNA editing enzyme catalytic subunit 3A; CXCL, C‐X‐C motif chemokine; CCL; C‐C motif chemokine; RSAD2, radical S‐adenosyl methionine domain‐containing protein 2; NEURL3, neuralized E3 ubiquitin‐protein ligase 3; USP18, ubiquitin‐specific peptidase 18; IFI44L, interferon‐induced protein 44‐like; CMPK2, mitochondrial UMP‐CMP kinase 2; LIPM, lipase member M; IFI, interferon‐alpha‐inducible protein; MX1, interferon‐induced GTP‐binding protein Mx1; IDO, indoleamine 2,3‐dioxygenase; MYEOV, myeloma‐overexpressed gene protein; ISG15, interferon‐stimulated gene 15; DDO, D‐aspartate oxidase; APO, apolipoprotein; GMPR, GMP reductase 1; DNASE2B, deoxyribonuclease‐2‐beta; ANKRD29, ankyrin repeat domain‐containing protein 29; INTU, protein inturned; SERPING1, plasma protease C1 inhibitor

FIGURE A4

Log2‐fold changes (logFC) of differentially expressed transcripts associated with M1‐shifted pro‐inflammatory responses and M2‐shifted anti‐inflammatory responses between ZIKV+ (P_ZV, NP_ZV, P_ILZV, and NP_ILZV) and ZIKV− (P_Ctr, NP_Ctr, P_IL3, and NP_IL3) PBMCs. CXCL, C‐X‐C motif chemokine; CCL; C‐C motif chemokine; IDO1, indoleamine 2,3‐dioxygenase 1; NR1H3, oxysterols receptor LXR‐alpha; CH13L1, chitinase‐3‐like protein 1; GM2, protein‐glutamine gamma‐glutamyltransferase 2; FCGR1A, high affinity immunoglobulin gamma Fc receptor I; EPAS1, endothelial PAS domain‐containing protein 1; IL, interleukin; TLR, toll‐like receptor; P2RY2, P2Y purinoceptor 2; SOCS, suppressor of cytokine signaling; STAT1, signal transducer and activator of transcription 1/alpha–beta; CD86, T‐lymphocyte activation antigen CD86; IRF5, interferon regulatory factor 5; BTK, Bruton's tyrosine protein kinase; HIF1A, hypoxia‐inducible factor 1‐alpha; NFKB1, nuclear factor NF‐kappa‐B p105 subunit; TNF, tumor necrosis factor; TGFB1, transforming growth factor beta‐1 proprotein; CD200R1, cell surface glycoprotein CD200 receptor 1; ARG1, arginase‐1

Differentially expressed immune‐related transcripts in ZIKV‐infected PBMCs from pregnant and non‐pregnant women

A total of 1965 genes were differentially expressed between ZIKV‐infected PBMCs from pregnant (P_ZV and P_ILZV) and non‐pregnant women (NP_ZV and NP_ILZV). Of these, 1001 genes (51%) were upregulated while 964 genes (49%) were downregulated in ZIKV‐infected pregnant conditions compared to non‐pregnant conditions. Seventeen and 21 of the 35 most overexpressed transcripts in ZIKV‐infected pregnant and non‐pregnant conditions, respectively, were protein‐coding genes. Higher logFC values of >1.5 for M1 marker genes, especially IFNγ (logFC = 2.1) were observed in ZIKV+ pregnant conditions. Conversely, higher logFC values for M2 marker genes, especially the macrophage mannose receptor 1 MRC1 (logFC = −6.6) and CD163 (logFC = −1.7) encoding the scavenger receptor cysteine‐rich type 1 protein M130, were observed in ZIKV+ non‐pregnant conditions (Figure 3).

FIGURE 3

Log2‐fold changes (logFC) of differentially expressed transcripts associated with M1‐shifted pro‐inflammatory responses and M2‐shifted anti‐inflammatory responses between ZIKV+ pregnant and ZIKV+ non‐pregnant conditions. BTK, Bruton's tyrosine kinase; CCL22, C‐C motif chemokine ligand 22; CD, leukocyte differentiation antigen (cluster of differentiation); CLEC10A, C‐type lectin domain‐containing 10A; HLA, human leukocyte antigen; IFNG, interferon gamma; IL, interleukin; IL1R, interleukin‐1 receptor; IRF4, interferon regulatory factor 4; KLF4, Krueppel‐like factor 4; MRC1, mannose receptor C‐type 1; NR1H3, oxysterols receptor LXR‐alpha; SOCS3, suppressor of cytokine signaling 3; STAT1, signal transducer and activator of transcription 1; TGM2, transglutaminase 2; TLR, toll‐like receptor; VEGFA, vascular endothelial growth factor A

Consistent with M1 marker gene observations, highly upregulated immune‐related transcripts in ZIKV+ pregnant relative to non‐pregnant conditions were pro‐inflammatory, including paraneoplastic Ma antigen family and allograft inflammatory factor 1 (PNMAL1; logFC = 8.2; Figure A5), mucin‐related transcripts MUC16 (logFC = 2.9; Figure A5) VWDE (logFC = 2.7), C‐C motif chemokine 4‐like CCL4L1 (logFC = 2.6), and IFNγ (logFC = 2.1; Figure 3). Consistent with IFNγ over‐expression, interferon‐induced responses were elevated in ZIKV+ pregnant relative to ZIKV+ non‐pregnant conditions. ZIKV+ pregnant conditions overexpressed transcripts for interferon‐induced proteins (IFI6, IFI16, IFI27, IFI35, IFI44L, IFIH1, IFIT1, IFIT3, and IFIT5; average logFC = 1.0 ± 0.3), interferon‐induced transmembrane proteins (IFITM1, IFITM2 and IFITM3; average logFC = 1.0 ± 0.5), and IFNγ‐activated signal transducer and activator of transcription proteins STAT1 (logFC = 0.6; Figure 3) and STAT3 (logFC = 0.3). The Ras‐related protein‐coding transcript RAB17 (logFC = 3.3) involved in the mitogen‐activated protein kinase pathway downstream to interferon activation was also among the most upregulated in ZIKV+ pregnant conditions (Figure A5). Only a few transcripts participating in interferon‐activated pathways, such as Janus kinase and microtubule‐interacting protein 2 JAKMIP2 (logFC = −0.8), as well as Ras guanyl‐releasing proteins RASGRP2 (logFC = −0.4) and RASGRP4 (logFC = −1.0), were downregulated in ZIKV+ pregnant compared to ZIKV+ non‐pregnant conditions.

FIGURE A5

Log2‐fold changes (logFC) of the 35 most differentially expressed transcripts between ZIKV+ PBMCs from pregnant (P_ZV and P_ILZV) and non‐pregnant women (NP_ZV and NP_ILZV). Solid bars represent regular gene transcripts, cross‐hatched bars represent pseudogenes and stippled bars long non‐coding RNAs. Abbreviations for regular gene transcripts: HLA‐DRB5, DR beta 5 chain of the HLA class II histocompatibility antigen; BIRC7, baculoviral IAP repeat‐containing protein 7; DNAAF3, dynein axonemal assembly factor 3; ZFP57, zinc finger protein 57; MRC1, macrophage mannose receptor 1; ADARB2, double‐stranded RNA‐specific editase B2; RP11‐497H16.6, β‐clucoronidase‐like protein SMA5; KRT73, type II cytoskeletal 73 keratin; RNF17, ring finger protein 17; AJAP1, adherens junction‐associated protein 1; VSIG4, V‐set and immunoglobulin domain‐containing protein 4; GNG8, guanine nucleotide‐binding protein G(I)/G(S)/G(O) subunit gamma‐8; SLC47A1, multidrug and toxin extrusion protein 1; MARCO, macrophage receptor MARCO; CH25H, cholesterol 25‐hydroxylase; OVCH1, ovochymase 1, KCNMA1, calcium‐activated potassium channel subunit alpha‐1; CRISP3, cysteine‐rich secretory protein 3; CEACAM3, carcinoembryonic antigen‐related cell adhesion molecule 1; MLXIPL, MLX interacting protein‐like; PNMAL1, paraneoplastic antigen‐like protein 8A; FAM203B, protein HGH1 homolog; SNX7, sorting nexin‐7; FAM21B, WASH complex subunit 2A; KRT1, type II cytoskeletal 1 keratin; BHLHA15, class A basic helix–loop–helix protein 15; AIF1L, allograft inflammatory factor 1‐like; TRBV10‐1, T‐cell receptor beta variable 10‐1; RAB17, Ras‐related protein Rab‐17; NEK2, serine/threonine‐protein kinase Nek2; HSPB7, heat shock protein beta‐7; RCAN2, calcipressin‐2; MUC16, mucin‐16; HBB, hemoglobin subunit beta; IGHA2, immunoglobulin heavy constant alpha 2; UTS2, urotensin‐2

The most downregulated transcripts in ZIKV+ pregnant relative to non‐pregnant conditions were associated with M1‐/M2‐type (HLA‐DRB5; logFC = −8.8) and M2‐type (MRC1, logFC = −6.6) macrophage responses (Figure A5 and Figure 3) (Ka et al., 2014). Another highly downregulated transcript, CEACAM1 (logFC = −2.6), encoded the carcinoembryonic antigen‐related cell adhesion molecule 1 that negatively regulates various immune pathways, including T‐cell proliferation, T‐cell responses, and cytokine production (Hosomi et al., 2013). Transcripts for receptors that facilitate viral entry into host cells, such as the ZIKV‐interacting lipoprotein receptor‐related protein 1 (LRP1; logFC = −0.8) and C‐type lectin domain family 4 member G (CLEC4G; logFC = −2.2) (Zhang et al., 2017), were also downregulated in ZIKV+ pregnant conditions.

Besides interferon pathways, autophagy and apoptosis pathways were also overexpressed in ZIKV+ pregnant conditions relative to ZIKV+ non‐pregnant conditions. Upregulated pro‐apoptotic transcripts included those encoding the apoptosis facilitator Bcl2‐like protein 14 (BCL2L14; logFC = 2.6), the apoptosis regulator Bcl2 (logFC = 0.4), Bcl‐2‐like protein 11 (logFC = 0.3), caspases 1, 3, 4, 7, 8, 10 (average logFC = 0.4 ± 0.1), and the Fas cell death receptor (logFC = 0.7). Another BIRC7 gene encoding pro‐/anti‐apoptotic baculoviral IAP repeat‐containing protein 7 was also among the most downregulated in ZIKV+ pregnant compared to ZIKV+ non‐pregnant conditions (Figure A5). Other downregulated transcripts included the PRAME gene encoding an anti‐apoptotic melanoma antigen (FC = −2.3) and the regulatory complex protein‐coding transcript LAMTOR4 (FC = −0.3) which inhibits autophagy by activating the mammalian target of rapamycin.

Differentially expressed osteoclast differentiation‐related transcripts across culture conditions

In addition to immune‐related transcripts, genes involved in osteoclast differentiation (KEGG ID: hsa04380) were also commonly overexpressed in ZIKV+ relative to ZIKV− PBMCs from either pregnant (P_ZV and P_ILZV compared to P_Ctr) or non‐pregnant conditions (NP_ZV and NP_ILZV compared to NP_Ctr). In pregnant conditions, 76 and 11 osteoclast differentiation genes were upregulated and downregulated, respectively. In non‐pregnant conditions, 67 and 12 osteoclast differentiation genes were upregulated and downregulated, respectively. The colony‐stimulating factor 1 (CSF1) gene essential for osteoclastogenesis was overexpressed in ZIKV+ compared to ZIKV− PBMCs in both pregnant and non‐pregnant conditions, while another essential gene encoding for the tumor necrosis factor‐related cytokine RANKL/TNFSF11 was overexpressed in ZIKV+ compared to ZIKV− PBMCs in pregnant conditions only (Boyle et al., 2003). Genes characteristic of the osteoclast lineage, including tartrate‐resistant acid phosphatase (ACP5), cathepsin K (CTSK), and β3 integrin (ITGB3), but not calcitonin receptor (CALCR), were upregulated in ZIKV+ relative to ZIKV− PBMCs from pregnant and non‐pregnant women. Differential expression analysis between ZIKV‐infected PBMCs from pregnant (P_ZV and P_ILZV) and non‐pregnant women (NP_ZV and NP_ILZV) showed the downregulation of 21 osteoclast differentiation genes and the upregulation of only eight genes in ZIKV+pregnant conditions relative to ZIKV+non‐pregnant conditions (Figure A6). However, many of these genes also participate in immune responses. This set of differentially expressed transcripts also did not include those essential for osteoclast differentiation, nor genes that typify the osteoclast lineage (except for CTSK).

FIGURE A6

Log2‐fold changes (logFC) of differentially expressed (a) osteoclast differentiation and (b) cardiomyopathy genes between baseline ZIKV+ (P_ZV and NP_ZV) and ZIKV− (P_Ctr and NP_Ctr) PBMCs. IFNG, interferon gamma; LILR, leukocyte immunoglobulin‐like receptor; STAT1, signal transducer and activator of transcription 1/alpha–beta; SOCS3, suppressor of cytokine signaling 3; AKT3, RAC‐gamma serine/threonine‐protein kinase; BTK, Bruton's tyrosine protein kinase; NFATC2, nuclear factor of activated T cells, cytoplasmic 2; CAMK4, calcium/calmodulin‐dependent protein kinase type IV; MAPK14, mitogen‐activated protein kinase 14; TGFBR, transforming growth factor beta receptor; SPI1, transcription factor PU.1; TYROBP, TYRO protein tyrosine kinase‐binding protein; SIRPA, Signal‐regulatory protein alpha; IFNGR1, interferon gamma receptor 1; NCF2, neutrophil cytosol factor 2; SIRPG/SIRPB1, signal‐regulatory protein gamma/beta‐1; OSCAR, osteoclast‐associated immunoglobulin‐like receptor; FCGR2A, low affinity immunoglobulin gamma Fc region receptor II‐a; FOS, proto‐oncogene c‐Fos; IL1R1, interleukin‐1 receptor type 1; CTSK, cathepsin K; TREM2, triggering receptor expressed on myeloid cells 2; ADCY5, adenylate cyclase type 5; ITG, integrin; DMD, dystrophin; TPM2, tropomyosin beta chain; PRKAG, 5′‐AMP‐activated protein kinase subunit gamma; ACE, angiotensin‐converting enzyme 2

In our separate analysis comparing gene expression between grouped infected (P_ZV, NP_ZV, P_ILZV, and NP_ILZV) and non‐infected (P_Ctr, NP_Ctr, P_IL3, and NP_IL3) cultures, osteoclast differentiation marker genes were similarly not differentially expressed, except for ITGB3 (logFC = 2.6). This suggests the potential effects of IL‐3 stimulation on the expression of osteoclast differentiation transcripts. We next investigated this by performing two differential expression analyses comparing non‐infected IL‐3+ and IL‐3− PBMCs from pregnant (P_IL3 vs. P_Ctr) and non‐pregnant conditions (NP_IL3 vs. NP_Ctr). 45 and 67 osteoclast differentiation‐related transcripts were upregulated in the P_IL3 and P_Ctr comparison and NP_IL3 and NP_Ctr comparison, respectively (Figure 4). Marker genes of osteoclast differentiation (CSF1, TNFSF11, ACP5, and ITGB3, but not CTSK) were among the most overexpressed in both P_IL3 and P_Ctr, as well as NP_IL3 and NP_Ctr comparisons (Figure 4).

FIGURE 4

Log2‐fold changes (logFC) of differentially expressed osteoclast differentiation transcripts between non‐infected IL‐3+ and IL‐3− PBMCs from (a) pregnant and (b) non‐pregnant women. Marker genes of interest are highlighted in blue. Abbreviations of marker genes: ACP5, tartrate‐resistant acid phosphatase; CSF1; colony‐stimulating factor 1; ITGB3, integrin subunit β3; TNFSF11, tumor necrosis factor superfamily 11

Differentially expressed cardiomyopathy‐related transcripts across culture conditions

Dilated cardiomyopathy (KEGG ID: hsa05414) and hypertrophic cardiomyopathy (KEGG ID: hsa05414) pathways were also commonly overexpressed in ZIKV+ relative to ZIKV− PBMCs from pregnant and non‐pregnant women (Figure 5). In pregnant conditions, 36 and 10 cardiomyopathy genes were upregulated and downregulated, respectively (Figure 5). In non‐pregnant conditions, 24 and 14 cardiomyopathy genes were upregulated and downregulated, respectively (Figure 5). Commonly upregulated transcripts encode for structural proteins, including cytoplasmic actin 1 and 2 (ACTB and ACTG1), tropomyosin alpha‐3 and ‐4 chains (TPM3 and TPM4), and prelamin‐A/C (LMNA; Figure 5). Other commonly overexpressed genes encode for integrins (ITGA5‐ITGA7, ITGAV, ITGB3, ITGB5, and ITGB8), adenylate cyclase type 3 (ADCY3), angiotensin‐converting enzyme (ACE), AMP‐activated protein kinase subunits gamma‐1 and ‐3 (PRKAG1 and PRKACA), cardiac‐type myosin‐binding protein C (MYBPC3), sodium/calcium exchanger 1 (SLC8A1), the endothelium‐derived vasoconstrictor peptides endothelin‐1 (EDN1), sarcoplasmic/endoplasmic reticulum calcium ATPase 2 (ATP2A2), voltage‐dependent calcium channel gamma‐8 subunit (CACNG8), transforming growth factor beta‐1 proprotein (TGFB1; Figure 5).

FIGURE 5

Log2‐fold changes (logFC) of differentially expressed cardiomyopathy‐related transcripts between ZIKV+ and ZIKV− PBMCs from (a) pregnant and (b) non‐pregnant women. MYH, myosin; ACT, actin; TNN, troponin; ITG, integrin; SLC8A1, sodium/calcium exchanger 1; EDN1, endothelin‐1; PRKAG, 5′‐AMP‐activated protein kinase subunit gamma; MYBPC3, cardiac‐type myosin‐binding protein C3; CACN, voltage‐dependent calcium channel subunit; LAMA2; laminin subunit alpha‐2; LMNA, lamin A/C; PLN, cardiac phospholamban; TPM, tropomyosin; ACE, angiotensin‐converting enzyme; DMD, dystrophin; PRKAC, cAMP‐dependent protein kinase catalytic subunit; TNF; tumor necrosis factor; TGFB1, transforming growth factor beta‐1; GNAS; guanine nucleotide‐binding protein G(s) subunit alpha isoforms short; ADCY3, adenylate cyclase type 3; PRKAB2/PRKAA1, 5′‐AMP‐activated protein kinase subunit beta‐2/alpha‐1; TTN, titin; ATP2A, sarcoplasmic/endoplasmic reticulum calcium ATPase; IL6, interleukin‐6

Differential expression analysis between ZIKV‐infected PBMCs from pregnant (P_ZV and P_ILZV) and non‐pregnant women (NP_ZV and NP_ILZV) predicted only six upregulated cardiomyopathy‐related genes and three downregulated genes in ZIKV+pregnant conditions relative to ZIKV+ non‐pregnant conditions (Figure A6). Nevertheless, other nervous and cardiovascular system‐related transcripts, including calcipressin‐2 (RCAN2; logFC = 3.1; Figure A5), UTS2 (urotensin‐2; logFC = 2.8; Figure A5), UTS2B (logFC = 2.3), and S100B (logFC = 2.6) were among the most upregulated in ZIKV+ pregnant conditions. On the other hand, muscle‐related transcripts including double‐stranded RNA‐specific editase B2 (ADARB2; logFC = −5.4; Figure A5), calcium‐activated potassium channel subunit alpha‐1 (KCNMA1; logFC = −2.6; Figure A5) and troponin T, slow skeletal muscle (TNNT1; logFC = −2.4) were among the most upregulated in ZIKV+ non‐pregnant conditions. Besides pregnancy status, in the absence of ZIKV infection, immune stimulation of PBMC cultures by IL‐3 alone in either pregnant (P_IL3 vs. P_Ctr) or non‐pregnant conditions (NP_IL3 vs. NP_Ctr) also led to the upregulation of eight cardiomyopathy‐related transcripts. These included ACTB, ACTG1, CACNG8, ITGA6, ITGB3, LMNA, SLC8A1, and TPM3.

DISCUSSION

In this study, we focused on PBMC responses to ZIKV to better understand general and specific immune responses to viral infections during pregnancy. We first predicted immune cell types from our PBMC transcriptomic data, then examined gene expression changes before and after ZIKV infection in pregnant and non‐pregnant conditions. We next identified differentially expressed genes between ZIKV‐infected pregnant conditions and ZIKV‐infected non‐pregnant conditions. Finally, we predicted differentially regulated non‐immune pathways across ZIKV infection, pregnancy, and IL‐3 stimulation conditions to identify candidate pathways that can affect pregnancy and birth outcomes.

Because PBMCs are heterogeneous and comprise different cell types that vary in gene expression (Whitney et al., 2003), our PBMC transcriptomic data enable a systems‐level perspective of gene expression profile changes between ZIKV‐infected and non‐infected, as well as non‐pregnant and pregnant conditions. Furthermore, to provide some cellular context to our data, we used the ABsolute Immune Signal (ABIS) deconvolution algorithm (Monaco et al., 2019) to predict proportions and mRNA abundances of immune cell types within PBMCs. Consistent with previous reports (Kieffer et al., 2017; Tilburgs et al., 2010), we found increased proportions and mRNA abundances of CD8 naïve T cells and CD4 memory T cells in uninfected pregnant relative to non‐pregnant conditions. This trend persisted even after ZIKV infection, and we speculate that baseline differential abundances in CD8 and CD4 T cells during pregnancy may affect immune responses against ZIKV and even pregnancy outcomes. ZIKV‐infected PBMCs showed lower abundances of CD4 memory T cells, and a statistically significant decrease in CD4 memory T‐cell counts was similarly reported in PBMCs infected with the Asian‐lineage ZIKV strain Nica 2–16 relative to uninfected controls (Michlmayr et al., 2017). Immune cell proportion and gene expression profile similarities between libraries NP_ILZV and NP_ZV and libraries P_ILZV and P_ZV indicated relatively smaller effects of IL‐3 stimulation on gene expression under ZIKV‐infected conditions. Predicted immune cell proportion data suggested that adding IL‐3 to ZIKV‐infected PBMCs promoted the differentiation of CD8 naïve T cells, but not the differentiation of CD4 naïve T cells. Similarly, a self‐reactive mechanism involving CD4+ IL‐3+ T cells was reported in an experimental autoimmune myocarditis mouse study (Anzai et al., 2019). In the study, IL‐3 produced by CD4 T cells promoted autoimmune inflammation, although IL‐3 is dispensable to CD4 T‐cell sensitization (Anzai et al., 2019). We speculate that a similar mechanism could occur during ZIKV infection and contribute to pro‐inflammatory responses, although these interactions will have to be confirmed by mechanistic studies. Compared to uninfected PBMCs, we observed lowered CD8 naïve T cells in ZIKV‐infected PBMCs relative to uninfected PBMCs. However, this trend was not supported by the Nica 2–16 study (Michlmayr et al., 2017). Previous studies also found monocytes to be the major component of ZIKV‐infected PBMCs (Foo et al., 2017; Michlmayr et al., 2017), but this was not captured in our study. Therefore, although our study highlights the utility of deconvolution algorithms in resolving immune cell populations in the absence of flow cytometry data, careful review and validation of the results are necessary. Furthermore, we were only able to make qualitative inferences from our deconvoluted data because of sample size limitations.

We performed differential expression analysis to identify gene expression changes before and after ZIKV infection in pregnant and non‐pregnant conditions. Compared to non‐infected controls, interferon‐inducible genes were commonly overexpressed in ZIKV‐infected PBMCs in both pregnant and non‐pregnant conditions. These genes overlapped with those previously reported in ZIKV‐infected PBMCs (Lanciotti et al., 2016), confirming our observations that the upregulation of interferon‐related genes during ZIKV infection was independent of IL‐3 stimulation. Pseudogenes and lncRNAs were also among the most up‐ and downregulated transcripts in ZIKV‐infected PBMCs compared to uninfected controls, as well as ZIKV‐infected PBMCs from pregnant conditions compared to non‐pregnant conditions. Because lncRNAs modulate the interferon antiviral response by affecting the activity of interferon‐stimulated genes and transcription factors, aberrant lncRNA expression has been observed in various viral infections (Qiu et al., 2018) and possibly even ZIKV infection from the same PRVABC59 strain used in our study (Ramaiah et al., 2016). The role of non‐coding genes, especially lncRNAs, in ZIKV infection thus warrants further study and can yield new insights on host–ZIKV interactions and pregnancy outcomes.

We detected greater upregulation of pro‐inflammatory M1 macrophage‐related transcripts and downregulation of anti‐inflammatory M2 macrophage‐related transcripts in PBMCs from pregnant women compared to non‐pregnant women 48 h post‐infection with the Asian‐lineage ZIKV strain Puerto Rico PRVABC59. Aggravated pro‐inflammatory monocyte responses during pregnancy were similarly observed with ZIKV strain MR766 (African lineage) infection (Foo et al., 2017) and influenza A virus infection (Le Gars et al., 2016). Nevertheless, our results contrast previous findings that PBMCs, particularly monocytes, infected with the Asian‐lineage ZIKV strain H/PF/2013 showed M2‐skewed immunosuppression responses (Foo et al., 2017). Hence, M1 vs M2‐skewed immune response may be strain‐specific and found in either African or Asian lineage. Because ZIKV can persist in serum for >60 days (Paz‐Bailey et al., 2019), the type of immune response may also depend on the infection stage during ZIKV isolation, similar to HIV‐1 coreceptor switch from M‐tropic to T‐tropic infection with disease progression (Connor et al., 1997). Further, PBMC infection rates ex vivo were found to vary across healthy donors (Messias et al., 2019), and thus, infection rate differences across strains, donors, or culturing conditions may result in different phenotypes. Because ABIS predicts only undifferentiated monocyte proportions, we were unable to directly compare M1/M2 marker gene expression with their corresponding macrophage abundances. There is no specified relationship between monocyte abundances and abundances of monocyte‐derived macrophages with M1‐like or M2‐like phenotypes (Italiani & Boraschi, 2014). Therefore, follow‐up studies with more detailed M1/M2 gene expression analysis would require further validation by flow cytometry and quantitative PCR approaches. Although previously reported pregnancy‐complication‐related fibronectin‐1 and ADAMTS9 genes (Foo et al., 2017) were not differentially expressed in our dataset, other upregulated transcripts in ZIKV+ pregnant conditions, including CX3CR1 (logFC = 0.4), IFITM, and even IFNγ, were associated with pregnancy defects in mice and humans (Buchrieser et al., 2019; Li et al., 2014; Saito et al., 1999). Likewise, pro‐apoptotic transcripts upregulated in ZIKV+ pregnant conditions may cause tissue damage and embryopathy (Mawson, 2016).

We also evaluated candidate non‐immune pathways that can potentially affect ZIKV‐associated pregnancy and birth outcomes. We observed the over‐expression of osteoclast differentiation genes with either IL‐3 stimulation or ZIKV infection. In line with gene expression similarities between libraries NP_ILZV and NP_ZV and libraries P_ILZV and P_ZV, our results suggest that ZIKV infection can induce IL‐3 activation. The role of IL‐3 on osteoclastogenesis is controversial, and both stimulatory (Hong et al., 2013) and inhibitory (Khapli et al., 2003) effects of IL‐3 on osteoclast differentiation have been reported. Nevertheless, based on our results, we were not able to conclude whether osteoclast formation was affected by pregnancy status, because differentially expressed genes in these pathways also participate in immune responses and are not marker genes characteristic of osteoclast differentiation or the osteoclast lineage. On the other hand, cardiomyopathy‐related genes appeared to be mainly stimulated by ZIKV infection. The acute phase of ZIKV infection has been associated with cardiovascular complications in humans, including mothers and infants born to ZIKV‐infected mothers (Minhas et al., 2017).

Collectively, our results reveal ZIKV infection‐associated and pregnancy‐associated differences in PBMC responses that could be of clinical and therapeutical relevance. As ZIKV and other infectious diseases continue to be a looming threat, larger‐scale integrative studies analyzing core and diverging immune responses against infection across microbial strains, as well as mechanistic studies disentangling molecular cross‐talk at the maternal–placental–fetal axis, will advance our current understanding on how microbial infections shape pregnancy outcomes.

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTIONS

Shen Jean Lim: Formal analysis (lead); software (lead); visualization (lead); writing‐original draft (lead); writing‐review & editing (supporting). Andreas Seyfang: Conceptualization (lead); funding acquisition (supporting); investigation (lead); methodology (lead); project administration (equal); resources (equal); supervision (equal); writing‐original draft (supporting); writing‐review & editing (equal). Samia Dutra: Investigation (supporting); project administration (supporting); resources (equal); writing‐original draft (supporting). Bradley Kane: Investigation (supporting); resources (equal); writing‐original draft (supporting). Maureen Groer: Conceptualization (lead); funding acquisition (lead); methodology (equal); project administration (equal); supervision (equal); writing‐review & editing (equal).

ETHICS STATEMENT

Research procedures adhered to the Declaration of Helsinki and were approved by the university Institutional Review Board (USF IRB protocol Pro00029803). Informed consent was obtained from all participants.