Differences in maternal gene expression in Cesarean section delivery compared with vaginal delivery

Cesarean section (CS) is recognized as being a shared environmental risk factor associated with chronic immune disease. A study of maternal gene expression changes between different delivery modes can add to our understanding of how CS contributes to disease patterns later in life. We evaluated the association of delivery mode with postpartum gene expression using a cross-sectional study of 324 mothers who delivered full-term (≥ 37 weeks) singletons. Of these, 181 mothers had a vaginal delivery and 143 had a CS delivery (60 with and 83 without labor). Antimicrobial peptides (AMP) were upregulated in vaginal delivery compared to CS with or without labor. Peptidase inhibitor 3 (PI3), a gene in the antimicrobial peptide pathway and known to be involved in antimicrobial and anti-inflammatory activities, showed a twofold increase in vaginal delivery compared to CS with or without labor (adjusted p-value 1.57 × 10–11 and 3.70 × 10–13, respectively). This study evaluates differences in gene expression by delivery mode and provides evidence of antimicrobial peptide upregulation in vaginal delivery compared to CS with or without labor. Further exploration is needed to determine if AMP upregulation provides protection against CS-associated diseases later in life.


Scientific Reports
| (2020) 10:17797 | https://doi.org/10.1038/s41598-020-74989-8 www.nature.com/scientificreports/ that CS delivery (with or without labor) compared with VD is associated with gene expression changes in the mother, which could potentially contribute to disease patterns later in life in the mother and offspring.

Results
Demographic and clinical data. A total of 324 mothers with RNA sequencing available were included in the analysis. Of these, 181 had a VD and 143 had a CS delivery (60 with and 83 without labor) ( Table 1). Mothers undergoing CS without labor were older than those undergoing VD or CS with labor (mean age CS without labor 34.3 years vs. CS with labor 31.4 years vs. VD 31.5 years; ANOVA p-value = 0.0002). Mothers undergoing CS without labor had a slightly lower gestational age at delivery than those having VD or CS with labor, although all deliveries were full term (mean gestational age CS without labor 38.8 weeks vs. CS with labor 39.1 weeks vs. VD 39.3 weeks; ANOVA p-value = 0.001). For mothers who had a CS, the indications for CS were also recorded ( Table 1). For CS with labor, the top reason for CS was failure to progress or prolonged labor (N = 34) and for mothers undergoing CS without labor, the top reason was repeat CS (N = 67).

CS (with or without labor) vs. vaginal delivery.
Samples from CS (with or without labor, N = 143) vs.
VD (N = 181) were compared. Ninety-one genes were found to be differentially expressed at FDR ≤ 0.05 and magnitude of log 2 (fold-change) ≥ 0.5 as highlighted in the volcano plot in Fig. 1A and listed in Supplementary  Table S1. Of the differentially expressed genes, 27 (30%) were upregulated in VD. Antimicrobial peptides were the top pathway and included PI3, which was the top gene upregulated in VD (Supplementary Table S2). No genes were found to be differentially expressed between delivery modes in an independent analysis of prenatal samples from these mothers, thus confirming the changes did not exist before delivery (lowest FDR > 0.99).     Table S4). Of these, 26 (63%) were upregulated in VD and 294 (37%) were downregulated. Table 3 lists the pathways obtained with at least 5 differentially expressed genes (separated by positive or negative fold change) in a pathway. Antimicrobial peptides and neutrophil degranulation were the top pathways for genes upregulated in VD. PI3 was the top gene upregulated in VD. No significant pathways were found for genes upregulated in CS without labor. No genes were found to be differentially expressed in an independent analysis of prenatal samples from these mothers (lowest FDR > 0.99).
CS with labor vs. without labor. Samples from CS with labor (N = 60) vs. CS without labor (N = 83) were compared. Differential gene expression analysis resulted in 263 genes at FDR ≤ 0.05 and log 2 (fold-change) ≥ 0.5 (Fig. 1D, Supplementary Table S5). Of these 263 genes, 60 (23%) were upregulated in CS without labor and 203 (77%) were upregulated in CS with labor. Table 4 lists immune pathways including interferon signaling and Cytokine signaling genes which were downregulated, and neutrophil degranulation and interleukin signaling genes, which were upregulated in CS with labor. No genes were found to be differentially expressed in the corresponding prenatal samples for the two groups (lowest FDR > 0.99).

Differentially expressed genes among all 3 groups. Likelihood Ratio Test (LRT) in the DESeq2 22
package was used for identification of genes that are not constant across the three delivery modes. The top 50 genes with the highest fold-change in either direction were selected from DEGs obtained for the two comparisons (VD vs. CS with labor, CS with labor vs. CS without labor) and all 41 DEGs were selected for CS without labor vs. VD. These genes were then intersected with 2,739 genes found to show change in expression across delivery modes with LRT (adjusted p-value < 0.001). This led to 81 genes (Supplementary Table S6) that were used for clustering samples across the three delivery modes (Fig. 2). PI3 was among the genes upregulated in VD compared to CS with or without labor (Fig. 3). While an independent analysis of prenatal samples did not show elevated expression in mothers undergoing VD (Fig. 3A,B), PI3 expression was increased in VD compared to CS in post-delivery samples (Fig. 3C,D).

Discussion
Differential gene expression was found between the delivery modes of VD, CS with labor, and CS without labor. Of particular interest, there is an enrichment of antimicrobial peptides (AMP) pathways in VD. PI3 gene, a known AMP, showed a twofold increase in VD (adjusted p-value < 10 -9 ) compared to CS. Genes in the AMP pathways are upregulated in VD samples compared to both CS with or without labor. AMPs are small molecular weight proteins that modulate host immune response against bacteria, viruses, and fungi 23,24 . Among the AMPs upregulated in VD compared to CS with or without labor (CTGS, DEFA1, DEFA1B, DEFA3, DEFA4, ELANE, GNLY, GZMH, LCN2, LTF, PI3, PRTN3, RNASE3, SEMG1), PI3 was the most significant in terms of adjusted p-value (< 10 -9 ) and fold change (> twofold). The peptidase inhibitor 3 (PI3) gene encodes elafin and demonstrates antimicrobial and anti-inflammatory activities 25 . Elafin has been shown to play a protective role in inflammatory bowel disease 26 , asthma 27 , and has been shown to be downregulated in acute stages of acute respiratory distress syndrome 28 . Elafin has also been proposed to be a modulator of innate immunity in the lower genital tract 29 . Furthermore, elafin expression was found to be elevated in the cervix as a response to pathogens during preterm labor 30 . Given the apparent protective antimicrobial and anti-inflammatory role of AMPs, we postulate that lower expression of PI3 in CS deliveries compared with vaginal deliveries, even if www.nature.com/scientificreports/ transient, may lead to may lead to an inflammatory and immune cascade and play a role in the reported immune and inflammatory adverse outcomes associated with CS delivery for both mothers and offspring 1-3,6-16 ; this warrants further exploration. Other differentially expressed pathways between delivery modes involve the inflammatory and immune pathways (neutrophil degranulation, interferon signaling). The significance of these pathways also require further investigation given the known stress and inflammation that occur both with labor 31 and with major surgery 17 .
Based on the known increased risk of chronic immune diseases in both mothers delivering by CS and their infants compared to VD, the finding of AMP upregulation, particularly PI3, in VD compared to CS has potential clinical implications. If further research confirms this association and a causal relationship is found, this could lead to the development of therapeutic interventions aimed to increase levels of AMPs in women undergoing CS.
Our main strength is that we included a large number of subjects in this cohort of pregnant women undergoing gene expression analyses from blood samples.
Because maternal blood sample collection was performed prenatally and postpartum, we were able to independently analyze samples from these two timepoints to conclude that changes in gene expression were influenced by the mode of delivery as no genes were differentially expressed across delivery modes in the prenatal samples.
A limitation of the study is that maternal prenatal blood was collected at various stages of the first and second trimesters between 12 and 27 weeks' gestation (mean: 26.4 weeks, 95% confidence interval: 25.9-26.8). It is possible that confounding factors contributing to the reason a mother had a CS vs VD, including a pregnant woman's medical and obstetric status and fetal status may lead to changes in gene expression. Although many confounding factors associated with delivery mode were adjusted for in the analysis (maternal age, maternal BMI, gestational age and computed ancestry), other factors were used as exclusion criteria for the analysis (multiple births, preterm births), and that no prenatal differences were seen in gene expression between mothers who then went on to have CS vs. VD, it is possible other confounding factors may have played a role in differences seen. It is known that following the rupture of membranes, there is an increased chance of ascending infections which may lead to increased AMP responses; in future studies, it may be valuable to account for if there was rupture of membranes, and duration between rupture of membranes and delivery. In addition, the study is limited due to the lack of long-term follow up of subjects to assess the development of possible adverse health outcomes.
In conclusion, differentially expressed genes were found in mothers according to delivery mode with enrichment of antimicrobial peptide pathways in vaginal delivery compared to both Cesarean section with or without labor. Further exploration is needed to confirm the differential expression seen, determine if similar changes are seen in gene expression, or gene products, in the offspring, and identify longer term epigenetic changes including associations with long-term clinical outcomes in both mother and baby given the risk of diseases associated with Cesarean section.  Differential gene expression and pathway analysis. EdgeR 36 was used for differential gene expression and previously published guidelines were followed for the analysis 37 . Genes with low counts were filtered out following EdgeR recommendation where at least n samples should have > 10/L counts per million where L is the smallest library size and n is the number of samples in the smaller group. In summary, lowly expressed genes were filtered and library sizes were recalculated, followed by trimmed mean of M-values (TMM) normalization (calcNormFactors), dispersion estimation using the negative binomial (NB) model (estimateDisp), extension of the NB model with quasi-likelihood (QL) methods (glmQLFit), and finally, testing for differential gene expression (glmQLFTest). Whole genome sequencing (WGS) data were also available for these subjects; common single nucleotide variants from WGS were used to compute principal components for depicting ancestry. Age, pre-pregnancy BMI, RNA sequencing processing batch, gestational age (in weeks), duration between birth and sample collection, and the top 3 principal components derived from genomic sequencing were used as covariates in the analysis. In order to detect statistically significant changes across a range of magnitude, a cutoff of absolute log 2 (fold-change) ≥ 0.5 and false discovery rate (FDR) ≤ 0.05 were used for selecting differentially expressed genes 38 . R package ggplot2 was used to create volcano plots (Fig. 1), and density plots for PI3 expression across delivery modes (Fig. 3). Hierarchical clustering was performed and heatmap plots (Fig. 2) were generated using heatmap.3 R package. Pathway analysis was performed using Reactome 39 , a manually curated pathway database that contains 12,608 human reactions organized into 2282 pathways involving 11,053 proteins as of version 71. Genes were separated by positive and negative fold-change in a given delivery mode for each pair-wise differential expression analysis and provided as input to Reactome for pathway analysis.
Statistical analysis. Comparisons of demographics and outcomes data were performed between women according to mode of delivery. Association between ethnicity and delivery mode was tested with Fisher's Exact Test, as implemented in the R 40 function fisher.test(). One-way analysis of variance (ANOVA) was performed using R function aov() to test for differences in age, BMI and gestational age across delivery modes. PI3 expression counts between delivery modes were compared using Student's t-test, as implemented in R function t.test().