Pregnancy is linked to faster epigenetic aging in young women

Data are derived from the Cebu Longitudinal Health and Nutrition Survey (CLHNS), a birth cohort study started in 1983 to 84 in Metropolitan Cebu, Philippines, and are available for download at https://dataverse.unc.edu/dataverse/cebu. Blank versions of relevant surveys are also provided in SI Appendix. The current study focuses on surveys conducted in 2005 (baseline) and 2009 to 2014 (follow-up). Both men and women in this study came from the original birth cohort, and so were born in the same geographical areas (barangays or neighborhoods) at the same time (during enrollment in 1983 to 84), but were not necessarily affiliated in any other way. Descriptive statistics of both baseline and follow-up samples are provided in Table 1. Surveys were administered in participants’ homes by a trained interviewer, and included questions about mental and physical health, behavior, sociodemographic context. A subset of questions focused on reproductive history, including number of known pregnancies, their duration, and outcomes. For men, the number of pregnancies fathered was self-reported. For women, pregnancy status in 2005 was reported at the time of sampling, and through back-calculation from subsequent surveys based on parturition within 9 mo of the original interview (all past pregnancies and their outcomes are recorded as part of the ongoing tracking process). All past pregnancies reported by participants were counted regardless of outcome, including miscarriages, stillbirths, and live births. Three women reported having had twins, which counted as a single pregnancy. Household income, parental education, and assets were used to create a composite score of socioeconomic status, described in more detail in ref. 30 and below. Follow-up samples focused on a subset of participants enrolled in additional surveys tracking new pregnancies between 2009 and 2014. Informed consent was obtained from all participants, and data collection was conducted with approval and in accordance with the Institutional Review Boards of the University of North Carolina at Chapel Hill, the Office of Population Studies Foundation, and Northwestern University.

A composite score of socioeconomic status was measured as a combination of income, education, and assets. Participants reported their annual income from all sources, including in-kind services, and the sale of livestock or other products by household members during the prior year, which were summed to determine total household income. Incomes were log-transformed. Maternal education (in years) was also reported. Participants also reported on ten assets (coded as 1 for possessing the asset in question, or 0 if not possessing the asset) that were selected to capture population-relevant aspects of social class, including electricity, refrigerators, air conditioners, color televisions, cable tv, tape recorders, electric fans, jeepneys, cars, trucks, and owning their residence. In addition, house construction type (i.e., light, mixed, permanent structure) was coded as 1, 2, or 3, respectively. Thus, asset scores ranged from 0 to 13. A principal components analysis was run on log income and assets, along with maternal education, at sample collection. The first dimension explained 70% of the variation in our composite SES score, and individual scores for the top component of variation were used as our measure of SES.

Urbanicity is associated with a range of health determinants in the Philippines (31). We used the urbanicity score developed by Dahly and Adair in this population to model this environmental effect (32). In brief, this scale was constructed using community-level data and validated with metrics of health in the CLHNS (32).

Principal components (PCs) of genome-wide genetic variation were included to control for potential population genetic structure. The derivation of these principal components has been described previously (33, 34).

Blood samples for DNA methylation were collected concurrent with in-home interviews. Baseline blood samples were collected in EDTA-coated vacutainer tubes from overnight fasted subjects. Follow-up blood samples were collected using capillary whole blood collected on filter paper. DNA was extracted using a standard protocol; 750 ng of genomic DNA was treated with sodium bisulfite (Zyme EZDNA, Zymo Research, Irvine, CA) and 160 ng of converted DNA was applied to the Illumina Infinium MethylationEPIC BeadChip under standard conditions (Illumina Inc., San Diego, CA). Technicians were blind to information regarding participant characteristics, and samples were randomly assigned to plate, chip, and row. Background subtraction and color correction were performed using Illumina Genome Studio with default parameters. Data were then exported into R for further analysis. Quality control for baseline was performed as part of a larger sample to confirm participant sex and replicate status. This was followed by quantile normalization on all probes including SNP-associated and XY multiple binding probes. To maximize the number of sites available for the epigenetic age calculator, probes with detection P-values above 0.01 were called NA for poor-performing samples only and were otherwise retained (25). The same quality control steps were followed for follow-up samples. Epigenetic age for Horvath, Hannum, PhenoAge, GrimAge, and DNAmTL clocks were calculated using the original online calculator (http://labs.genetics.ucla.edu/horvath/dnamage/–new users are directed to https://dnamage.clockfoundation.org/). Principal component-based epigenetic clocks were calculated using the methods and code described in ref. 35. Background-corrected DNAm beta values (which represents the ratio of methylated beads to methylated and unmethylated beads) were processed further using the calculator’s internal normalization algorithms. The DunedinPACE clock was generated using the DunedinPACE calculator, available at https://github.com/danbelsky/DunedinPACE. Clock values reported in Table 1 are residualized on chronological age and so are centered on zero (except DunedinPACE, which represents estimated pace of aging in units of years). Equivalent months of epigenetic aging was determined by multiplying the SD of each epigenetic aging measure in the sample by the standardized effect size for that measure. For aging measures with units in years, this value was multiplied by 12 to provide equivalent months of epigenetic aging.

We first examined cross-sectional DNAm measures of aging for women at baseline who had never been pregnant compared to those who had. We next examined the relationship between cross-sectional DNAm measures of aging and gravidity as a continuous variable. In both cases, epigenetic age was the outcome of interest, with age, gravidity, pregnancy status at the time of blood sampling (i.e., whether or not a woman reported being pregnant at the time of the blood sample or was retrospectively deemed to be pregnant based on subsequent follow-up surveys), composite score of socioeconomic status, an urbanicity score (32), smoking status, and the top-10 principal components of genetic variation as covariates. Because the epigenetic clocks that we examined all aim to quantify potentially overlapping dimensions of biological aging, we did not consider tests of each clock to be independent tests and did not correct for multiple comparisons. This approach was prespecified in our preregistered analysis plan (https://osf.io/mqb37).

Exploratory data analysis and quality checks revealed potential outliers with values greater than three SD from the mean for several epigenetic clocks (SI Appendix, Fig. S1). In accordance with our preregistered analysis plan (https://osf.io/mqb37), we first fit ordinary least squares models including these values and examined diagnostics plots. These extreme observations were often high-leverage with large Cook’s distances, and had disproportionate influence on model estimates. The mean number of pregnancies reported were similar for those with epigenetic age estimates >3 SD from the mean compared to those with epigenetic age estimates =< 3 SD from the mean (mean pregnancies 0.57, SD 0.85 vs. mean pregnancies 0.52, SD 0.87, respectively). However, to moderate the influence of any single data point on the final estimates, and again in accordance with our preregistration protocol, we fit all models using a robust regression method defined by Yohai (36). We base our conclusions on the outcomes of the robust regression models. However, we also provide estimates from the ordinary least squares regression (including outliers) in SI Appendix, Tables S1–S20.

Next, we modeled longitudinal effects of reproduction by examining whether changes in DNAm measures of aging were associated with changes in gravidity. Here, change in epigenetic age was the outcome of interest, change in gravidity was the predictor of interest, and change in age, baseline pregnancy status at the time of blood sampling, composite score of socioeconomic status, urbanicity score, smoking status, and the top-10 principal components of genetic variation were included as covariates. Epigenetic clock estimates are based on DNAm derived from blood leukocytes. Because DNAm plays a key role in cellular identity, immune cell composition at the time of sampling may partly confound estimates of epigenetic age. To address this, we ran sensitivity analyses on all models that included bioinformatically estimated proportions of CD4T, CD8T, natural killer (NK), B cell, monocytes, and granulocytes (37). Cross-sectional analyses at baseline included baseline measures of immune cell counts, while longitudinal models included changes in cell counts between baseline and follow-up. Cross-sectional models using the same covariates (with the exclusion of current pregnancy status), including sensitivity models incorporating bioinformatically estimated cell composition, were run in same-aged cohort men. In keeping with our preregistration plan, all models were also fit using principal component-based versions of Horvath, Hannum, PhenoAge, GrimAge, and DNAmTL clocks (35). Data for the cross-sectional analysis of women is provided in SI Appendix, Data File 1. Data for the longitudinal analysis of women is provided in SI Appendix, Data File 2. Data for cross-sectional analysis of men is provided in SI Appendix, Data File 3.

We first tested whether women who had been pregnant by early adulthood appeared biologically older than women who had not been pregnant. Following our preregistered analysis plan (https://osf.io/mqb37), we fit robust models as defined by Yohai (36) with pregnancy (ever vs. never pregnant) as the exposure and six epigenetic clocks as the outcomes of interest. Epigenetic clocks included the Horvath and Hannum first-generation clocks trained on chronological age; the PhenoAge and GrimAge second-generation clocks trained on mortality; the DunedinPACE pace of aging measure trained on the change in 19 indicators of organ-system integrity across two decades; and DNAmTL, a DNAm surrogate measure of leukocyte telomere length. As covariates, we included a composite measure of socioeconomic status that included all sources of household income, education, and family assets reflecting population-relevant aspects of social class; a measure of urbanicity of the participants’ primary residence; pregnancy status at the time of the blood sample; smoking status at the time of blood sample; and the top ten principal components of genome-wide genetic variation.

For all DNAm measures of aging, women with a history of at least one pregnancy appeared biologically older than women who had never been pregnant (Table 2). The effect sizes ranged from 0.14 SD for the GrimAge clock to 0.28 for DunedinPACE (Table 2). In contrast with other epigenetic clocks, pregnancy was associated with shorter DNAmTL, consistent with more advanced biological age. The effect of having been pregnant was equivalent to between 4.0 and 14.2 mo of accelerated aging for first- and second-generation clocks, an accelerated pace of aging of 3.3% per year according to DunedinPACE, and a shortening of 0.03 kilobases according to DNAmTL. Results for both robust models and ordinary least squares models including high-leverage data points are provided in SI Appendix, Table S1.

Next, we tested whether women who experienced a greater number of pregnancies appeared biologically older than women with fewer or no pregnancies. Following our preregistered analysis plan (https://osf.io/mqb37), we fit robust models as defined by Yohai (36) with gravidity (number of pregnancies) as the exposure of interest, and the same six DNAm epigenetic clocks as the outcomes of interest. Socioeconomic status, urbanicity, pregnancy status, genetic variation, and smoking status were again included as covariates. For all DNAm measures of aging, women with a history of more pregnancies appeared biologically older than women with fewer pregnancies (Table 3). The per-pregnancy effect sizes ranged from 0.07 SD for the GrimAge clock to 0.16 SD for DunedinPACE. These effects are equivalent to 2.0 and 5.2 mo per pregnancy for first- and second-generation clocks, an accelerated pace of aging of 1.9% per year per pregnancy according to DunedinPACE, and shortening of 0.011 kilobases per pregnancy according to DNAmTL. Results for both robust models and ordinary least squares models are provided in SI Appendix, Table S3. Gradients in biological aging by number of pregnancies at study baseline are shown in Fig. 1.

DNAm in blood can be affected by leukocyte composition at the time of sampling (38), which may be reflected in DNAm measures of aging (39, 40). To test the robustness of our findings to differences in leukocyte composition between women, we refit the above models with the addition of estimates of CD4T, CD8T, natural killer cells, B cells, monocytes, and granulocytes (37). Similar to models that did not include cell count, women with a history of at least one pregnancy appeared epigenetically older than women with no pregnancies, and women with a greater number of pregnancies appeared epigenetically older than women with fewer pregnancies. Standardized coefficients for the effect of gravidity on DNAm aging increased with the inclusion of cell count estimates for all clocks except Horvath and DNAmTL, which decreased modestly. Results for both robust models and ordinary least squares models controlling for cell counts are provided in SI Appendix, Tables S2 and S4. Cross-sectional models were also run using principal component-based versions of all clocks and yielded similar results (SI Appendix, Tables S11–S14).

We next asked whether the pace of biological aging that women experienced during the follow-up period was related to the number of pregnancies a woman had experienced during that period. By comparing women at follow-up to themselves at baseline, this longitudinal approach minimizes potential confounding by time-invariant, between-individual variation in health, access to resources, and other social and environmental exposures (29). Women who had more pregnancies between baseline and follow-up showed greater changes in both Horvath and Hannum first-generation clocks, but no significant changes according to the other DNAm measures of aging (Table 4 and Fig. 2). Horvath and Hannum effect sizes were similar to those of our cross-sectional analysis. Horvath epigenetic age increased 0.06 SD—equivalent to 2.9 mo—for each additional pregnancy (95% CI b = 0.01 to 0.11, P = 0.024). Hannum’s epigenetic clock increased 0.06 SD—equivalent to 2.4 mo—for each additional pregnancy (95% CI b = 0.01 to 0.11, P = 0.026). Results for both robust models and ordinary least squares models are provided in SI Appendix, Table S5. Sensitivity analyses that included changes in immune cell composition yielded similar results, with standardized coefficients showing almost no change with the inclusion of cell count estimates. Results for both robust models and ordinary least squares models controlling for cell counts are provided in SI Appendix, Table S6. Longitudinal models were also run using principal component-based versions of all clocks and yielded similar results (SI Appendix, Tables S15 and S16).

In our final analysis, we tested whether the number of pregnancies cohort males fathered was associated with accelerated biological aging (n = 910). The men in this sample were born in the same geographical areas at the same time, allowing us to further disentangle potential sociocultural confounders of reproductive activity from the direct energetic investment of pregnancy and breastfeeding in the CLHNS. Using the same cross-sectional models described for the women in this cohort (minus the covariate for current pregnancy status), we found no evidence that men who reported fathering at least one pregnancy were biologically older than men who did not report fathering any pregnancies (SI Appendix, Table S7). In fact, men who fathered at least one pregnancy looked biologically younger according to Levine’s PhenoAge clock than men who did not father any pregnancies, but this effect was not statistically significant in sensitivity analyses that included immune cell composition (SI Appendix, Table S8). When examining number of pregnancies as a continuous variable, we found no evidence that men who reported fathering a greater number of pregnancies were epigenetically older than men who fathered fewer pregnancies (SI Appendix, Table S9). Sensitivity analyses that included immune cell composition yielded similar results, with standardized coefficients showing little change with the inclusion of cell count estimates (SI Appendix, Table S10). Cross-sectional models for men were also run using principal component-based versions of all clocks and yielded similar results (SI Appendix, Tables S17–S20).

C.P.R., N.R.L., D.B.C., and C.W.K. conceptualized and designed the research; C.P.R., N.R.L., D.B.C., J.L.M., D.T.S.L., P.A., and C.W.K. collected, curated, and processed the data; M.K. contributed reagents or analysis tools; C.P.R., D.W.B., and C.W.K. analyzed data; C.P.R., N.R.L., D.B.C., J.L.M., D.T.S.L., P.A., D.W.B., M.S.K., and C.W.K. read and approved the final manuscript; and C.P.R. and C.W.K. wrote the paper.

D.W.B. is listed as an inventor on a Duke University and University of Otago invention, DunedinPACE, that was licensed to a commercial entity.