Epigenetic aging and lifespan reflect reproductive history in the Finnish Twin Cohort

Domain: biology
Published: 2026-01-08
Editorial outcome: major_revision

Abstract
Reproductive history is closely linked to health, yet its relationship with biological aging and survival remains uncertain. We investigated this in the Finnish Twin Cohort, a population-based study that enables modeling of full childbearing history while controlling for common risk factors, through questionnaires and civil registries. We model the association between reproductive trajectories and survival in 14,836 women, and assess biological aging in a subset of 1054 participants using the PCGrimAge, an algorithm trained to predict biological aging and mortality risk from DNA methylation. We identify six distinct reproductive trajectories describing different timing and number of childbearing events. Women with the most live births throughout their lives (mean 6.8, SD 2.4) and nulliparous women showed accelerated aging and elevated mortality risk. These findings support the disposable soma theory of aging in modern humans, and provide valuable insights into the genetic and lifestyle-related determinants of lifespan.

Introduction
Life history theory suggests that, within a single lifespan, organisms must balance limited resources between reproduction and survival1. Consequently, aging — understood here as the gradual decline in physical capabilities over time — is viewed as a consequence of resource allocation favoring reproduction over bodily maintenance, as proposed by the disposable soma theory of aging2,3. Multiple studies have found a U-shaped relationship between reproduction and mortality in humans and others species: both nulliparity and giving birth to many offspring have been linked to increased mortality risk4–6. At one end of the life-history spectrum, lower survival prospects and low reproductive output may both be caused by poorer initial health - whereas at the other end of the spectrum, shorter lifespan in individuals with high numbers of offspring may reflect the physiological implications of a greater reproductive investment5. Interestingly, earlier reproduction has also been linked to a shorter lifespan, which may reflect a coupling between reproduction and survival independent from the number of offspring7. Yet, to our knowledge, most studies have only focused on the effects of single parameters, e.g. age at first reproduction, on lifespan5–8, and have not assessed the timing of each consecutive childbirth. But a reproductive history (including both the number and the timing of childbearing events) is highly multivariate: while classical modeling approaches are based on selected indicators (typically offspring count, age at first breeding or average interbirth interval), they need to balance parameter space, collinearity between parameters, and accuracy in the representation of reproductive trajectories. Mixture modeling methods, such as latent class analysis (LCA)9 classify reproductive histories into categories based on timing and frequency, capturing nuanced patterns while reducing overfitting and minimizing oversimplifications inherent in tabular methods. At the same time, lifespan is only a proxy for the aging process itself, and can be influenced by external events: as a consequence, it is still unknown whether increased reproductive output really accelerates the aging process itself, what is the effect of the timing of childbirth, and whether any age acceleration is already observable before old age.

DNA methylation-derived age estimates—the epigenetic clocks— can be used to approximate one’s biological age as the apparent age of the body10. Furthermore, the difference between known chronological age and estimated biological age gives us access to an individual’s aging rate, i.e. the relative speed of each individual’s own clock11. Epigenetic clocks can be categorized based on their training focus: those trained to estimate chronological age as closely as possible (such as Horvath’s clock12 and Hannum’s clock13), and those trained to assess phenotypic age (including the DunedinPACE14, PhenoAge15 and GrimAge16 clocks). The GrimAge epigenetic clock, in particular, excels in predicting time-to-death, offering a valuable measure of age-related physiological impairment, years before the end of life16. A recently developed principal component (PC)-based version of the GrimAge clock (hereafter PCGrimAge), uses PCs derived from CpG-level DNA methylation data for biological age estimation17. The PCGrimAge has enhanced the reliability of biological age estimates, as shown by the reduced deviations between replicates compared with the original CpG-level clock17. The present study is based primarily on that clock.

Some pioneering cross-sectional studies18–20 have linked the number of pregnancies to epigenetic age acceleration, but the evidence has been equivocal21. Recently, Ryan et al. (2024)18 found that increasing parity accelerated epigenetic aging in 22 year old women using six independent epigenetic clocks. However, like in other previous studies, the focus is restricted to younger women without assessing their survival prospects severely limits our understanding of the long-term consequences of reproduction on aging. What is more, to our knowledge, previous studies have focused on the number of children or the age at first reproduction, regardless of the timing of subsequent pregnancies, and thus could not consider the effects of the full childbearing history, including later-life pregnancies18–20. Finally, it has so far been challenging to assess the influence of confounding factors, such as early life, socio-economic background, or genetic determinants, that could affect both reproduction and lifespan: for example, Long and Zhang recently identified 98 variants that promoted reproduction but decreased survival7. Ryan et al. controlled for genetic background through a longitudinal study design, where they found that pregnancies between age 25 to 31 years accelerated epigenetic aging using two epigenetic clocks, yet whether these effects were transient or persisted later in life remains unknown18. Thus, the long-term associations of epigenetic aging trajectories and survival with lifelong reproductive history on the one hand, and familial and lifestyle factors on the other remain largely unknown.

The current study addresses this knowledge gap by investigating whether childbearing associates with accelerated aging and lifespan in modern day humans. Based on a dataset consisting of 14,931 Finnish twin women22, we assess the effect not only of the number of offspring, but also the timing of each childbirth. To that aim, we use a dimension-reduction approach based on latent class analysis (LCA)9: briefly, we identify K typical reproductive trajectories based on a very large reference sample (n=14,931) in our focal population (19-20th century Finland) and describe each individual as a weighted mixture of these K components. This allows for (1) the representation of any individual trajectory as a fixed-length vector of weights of size K, and (2) for analysis of the outcomes of each of the K typical strategies without attempting to separate the interlinked effects of e.g. age at first reproduction and total number of offsprings, and (3) to explore whether our study participants can be differentiated into subgroups based on their childbearing history, simultaneously providing valuable insights into the reproductive behavior of Finnish women born since the 1880’s until the mid-20th century. Given the potential balance in resource allocation between reproduction and somatic maintenance, we expect that women with higher lifetime reproduction trajectories will display accelerated epigenetic aging and a shorter lifespan. We test this by modeling survival as a function of the reproductive trajectories. Further, among women with methylation data (n=1054), we assess differences in biological age acceleration using the PCGrimAge16,17 algorithm. As we studied twin individuals sampled as twin pairs, we accounted for the non-independence between twins from the same family by including twin pair ID as a random effect in our models. Our method also incorporates data on socioeconomic background and lifestyle-related factors, adjusting for prevalent risk factors of later life morbidity of the 20th century, allowing a rigorous examination of the interplay between reproduction, aging and lifespan.

Results
Latent class analysis identifies six reproductive classes
In order to account not only for the number but also for the timing of births, we first modeled individual reproductive trajectories as a mixture of K typical trajectories inferred from lifelong live birth data (Methods, Supplementary Results) using latent class analysis (LCA9). This data included 10,783 women, who had had at least one live birth during their life. The distribution of ages at childbirth in the raw data is shown in Fig. 1. We tested models with live births summed over 1 to 4 years bins, and K=1 to 8 latent classes, and selected the one minimizing sample-size adjusted BIC, with the constraint of keeping a meaningful sample size in each class (at least 3% of women with maximum mixture weight in that class), and maximizing entropy and class-specific posterior probabilities (Supplementary Table 1). In the selected model, individual reproductive trajectories were described as a mixture of K=6 classes (hereafter reproductive classes), where live births were summed over 3-year age bins (Methods; latent class parameters are given in Table 1, cross-class classification probabilities are given in Supplementary Table 2). In this framework, each reproductive class is thus defined as a typical timing and number of live births (Fig. 2a): for example, a woman in class 1 had, on average, 2 birth givings (SD=1) and a high probability to give birth before age 24, while a woman in class 6 would have many live births (average 6.8, SD=2.4) interspersed throughout her reproductive life. This result was well supported by parametric bootstrapping (Supplementary Fig. 1); under this model, women are mostly assigned to a single class (Supplementary Fig. 2). Women who had not given birth (n=4148) were added to this data and assigned to a separate class (class 0, hereafter nulliparous) with a posterior probability of ~1 for further analyses.

Reproductive classes predict survival
Survival differed significantly between reproductive classes when modeled as mixed Cox proportional hazards, weighted on probabilities of belonging to each class (n=14,836, Methods, Fig. 2b, c). All of our survival models were adjusted for the historical birth cohort to control for differences in living circumstances in Finland during the range of birth years in our sample (Methods, Supplementary Table 3, Supplementary Results). When accounting for historical birth cohort and relatedness, the nulliparous women, or with few live births early in life (classes 0, 1 and to a lesser extent 2) had increased risk of death compared to the reference class i.e. class 3 with hazard ratios (HR) of 1.43 (z=7.08, p=4×10−12), 1.18 (z=2.78, p=5.3×10−3), and 1.10 (z=1.81, p=7×10−2), respectively, see 95% confidence intervals (CI) in Fig. 2b). Symmetrically, women with many live births throughout life (class 6 with 6.8 births on average, Table 1) also had increased mortality (HR=1.25, z=2.75, p=6×10−3, Fig. 2b). This pattern was attenuated but not abolished when correcting for known risk factors body mass index (BMI) tobacco and alcohol use, and education (Fig. 2c, Supplementary Table 3, Methods), also for the complete case sample with no covariate imputation (Supplementary Fig. 3, Methods). We used alternative modeling strategies to assess the robustness of our results. First, we implemented a singleton approach by drawing 100 random subsamples, each including one individual per twin pair from the 9496 families, and reran our survival models on these samples. Second, we analyzed survival time using fixed-effects models with cluster-robust standard errors. Both approaches consistently showed elevated hazard ratios for nulliparous women, early mothers, and the high lifetime reproduction group (classes 0, 1 and 6). The increased risks for nulliparous women and the high lifetime reproductive group remained statistically significant after adjusting for risk factors, underscoring the robustness of these associations (Fig. 3, Supplementary Tables 4 and 5, Methods: Alternative analysis strategies).

Reproductive classes predict aging
Similarly to survival, there were significant differences between the reproductive classes in epigenetic aging measured by the PCGrimAge clock in a subsample of 1054 women (Fig. 2d, e, see statistics on pairwise comparisons between all classes in Supplementary Table 6). In our minimally adjusted model (controlling for chronological age and relatedness, Methods), we observed the most accelerated aging for women with few live births in early life (class 1) compared with the women giving birth in their late 20s and early 30s (classes 3-5), with 1.42 years (SE=0.55, p=0.01) higher biological aging rate than class 4, that had the lowest aging rate, and 2.4 live births on average (see class-specific estimates and 95% CI in Fig. 2d). When additionally adjusting for common risk factors (BMI, tobacco and alcohol use, and education), the most striking difference in epigenetic age acceleration was observed between the class with the highest reproductive output (class 6) and the class 5, which was among the lowest reproductive investment (2.0 live births on average, and had the latest timing of reproduction): class 6 had 1.35 years (SE=0.54, p=0.01) higher epigenetic aging rate than class 5 (Fig. 2e). We also found accelerated aging for the nulliparous women, who had 1 year (SE=0.38, p=8×10−3) more accelerated biological aging compared to class 5, the latest reproduction group (Fig. 2e). These results were also supported by our singleton approach with 100 random draws of individuals from each of the 611 families in our epigenetic age acceleration analyses (Fig. 3, Supplementary Table 4). We tested whether we would observe similar patterns with two other published epigenetic clocks that are trained to capture different aspects of biological aging: DunedinPACE14, which is trained to predict pace of biological aging using 19 indicators of organ-system integrity, and PC-based17 PhenoAge, a clock that is developed to measure phenotypic age from clinical measures15. In both, the unadjusted and adjusted DunedinPACE models, we observed the highest pace of aging for the class with most birth givings (class 6), compared to the nulliparous and classes 2–5 (Supplementary Fig. 4, Supplementary Table 6). We also found accelerated aging for the early mothers (class 1), yet this difference did not remain significant in the risk-factor-adjusted models. In contrast to our results for PCGrimAge, we did not observe accelerated aging for nulliparous women using DunedinPACE, while analysis with PCPhenoAge showed no differences between the classes (Supplementary Fig. 4, Supplementary Table 6).

Parity, age at first child, and later health
Our latent class results are consistent with tabular approaches, as we modeled survival and epigenetic age acceleration as a function of parity and age at first childbirth using models with same covariates as in the fully adjusted models (Methods). We observed increased mortality risk for the nulliparous women (HR=1.37, z=5.7, p=1.2×10−8, 95% CI=1.23–1.52), and women with five or more live births (HR=1.22, z=2.78, p=5.4×10-3, 95% CI=1.06-1.40) compared to the women with three live births (used as a reference), yet no other differences between parity groups were observed (Fig. 4a). Also, the nulliparous women and women with five or more live births had accelerated aging by PCGrimAge and DunedinPACE, but only nulliparous were statistically different from other groups in the risk factor adjusted models (0.95 years of age acceleration for nulliparous compared to women with 1 live birth, SE=0.27, p=4.9×10−3, Fig. 4b, Supplementary Fig. 5). Later age at first childbirth was associated with slightly decreased mortality risk (HR=0.98, z=−4.12, p=3.7×10−5, 95% CI=0.97–0.99) and lower DunedinPACE of aging (Estimate=-1.5×10−3, SE=6.7×104, p=2.5×10−3) in our unadjusted models, yet in risk-factor-adjusted models this effect was statistically nonsignificant.

Methods
This study has been approved by the appropriate national research ethics committees, with the most recent approval granted by the Hospital District of Helsinki and Uusimaa ethics board in 2018 (#1799/2017). Blood samples for DNA analyses were collected from each participant after they had signed a written informed consent. The samples and data used for the research were obtained as part of the Finnish Twin Cohort studies and have been deposited in the Biobank of the Finnish Institute for Health and Welfare in accordance with the Finnish Biobank Act (2013). The authors state that all procedures involved in this study adhere to the ethical standards of human experimentation, as well as to the principles of the Declaration of Helsinki.

Software
All analyses were performed using Mplus (version 8.2 used throughout the manuscript)42 and R (version 4.2.3 used throughout the manuscript)43 softwares. Specific R packages used in the analysis are detailed in the corresponding method sections. All plots were generated in R with the package ggplot244 and further edited using Inkscape (version 1.4.2)45.

Participants
The study participants were a subset of the older Finnish Twin Cohort (FTC)22, a population-based prospective cohort consisting of twins from monozygotic (MZ) and dizygotic (DZ) pairs born before 1958, with both twins alive in 1967. In 1975, a mailed questionnaire survey was sent to all twin pairs with known residential address22. The questionnaire data includes self-reported information on health and health-related factors. This included information on family members, education, weight, height, and lifestyle-related factors such as smoking and alcohol consumption. The current study includes 14,931 women from the cohort, who were born between 1880 and 1957, and had a known number of live births. The full sample was used to compile the reproductive trajectories in latent class analysis, while a sample of 14,836 women who lived at least 40 years was used to assess the link between reproductive history and mortality risk. For a subsample of 1054 participants, blood samples were taken during 1994–2020. These samples were used to assess genome-wide DNA methylation and epigenetic age acceleration.

Live birth data
Live birth data was curated for all the study participants from the 1975 questionnaire. Participants were requested to report the birth years of each son and daughter. In case a woman had reported multiple childbirths within the same calendar year, these were considered to be from a multiple pregnancy (i.e. twins, triplets, etc.). This data was validated for all participants born between 1950–1957 (n=4431, 41.1%) by retrieving the exact birth dates of each child from the Digital and Population Data Services Agency on 31 December 2009.

Contextual data
We adjusted our statistical analyses with variables that were hypothesized to contribute to variance in epigenetic age acceleration and mortality risk. See sample sizes by variable categories in Supplementary Table 3.

Historical birth cohort. Our study participants were born between 1880 and 1957, a period during which Finland went through drastic societal and political changes, with likely influence on the study participants’ health and survival prospects, as well as reproductive opportunities. To account for these differences in living circumstances, we assigned each participant to their birth cohort, reflecting differing living circumstances in each time period: 1880–1898 Industrial era, 1899–1917 Pre-independence, 1918–1938 Post-independence, 1939–1945 World War II, 1946–1957 Cold War.

BMI. Body mass index (BMI, kg m−2) has been shown to be associated with epigenetic aging and mortality risk46,47. For survival analysis, BMI was calculated from self-reported weight and height in 1975, which have been previously shown to correlate well (r(264)=0.95, CI=0.94–0.96) with measured values for women of the same cohort38. For epigenetic aging analyses, BMI at the time of blood-draw was used. As BMI may have a non-monotonous relationship with health47, we used centered BMI as a linear term as well as a quadratic term in both survival and epigenetic aging models.

Smoking. Smoking is also a known risk factor for mortality and has been shown to influence DNA methylation and accelerate epigenetic aging46,48. For survival analysis, participants were attributed a smoking status (never, former, occasional (non-daily), light (1–9 cigarettes per day), medium (10–19 cigarettes per day), heavy 20 or more cigarettes per day) based on questionnaire data from 197549. At the time of blood sampling, the participants were asked a three-level smoking status (never, former, current), which was used as a covariate in the epigenetic aging analysis. We validated the reliability of self-reported smoking information by calculating correlations between self-reported smoking pack-years (calculated as packs smoked per day multiplied by years as a smoker), preprocessed AHRR DNA methylation at locus cg05575921, which has been previously found to be hypomethylated in heavy-smokers50, as well as blood-based cotinine measures51. The Pearson correlation coefficient between self-reported pack-years and preprocessed cg05575921 methylation was r(1660)=−0.48 (p<2.2×10−16), whereas the correlation between pack-years and plasma cotinine was r(224)=0.61 (p<2.2×10−16).

Alcohol consumption. Alcohol use is associated with DNA methylation, epigenetic aging, as well as morbidity and mortality risk52,53. An alcohol use index was compiled from the 1975 questionnaire, first as a continuous variable reflecting average daily use in grams of ethanol49, further binned into a categorical variable (abstainer, infrequent (<1.3g per day), low (<25g per day), medium (<45g per day), high (<65g per day), higher (≥65g per day) as further described in Zhao et al. (2023)54. As heavy drinking was rare among the women in our cohort, the medium, high, and higher categories were combined into a medium or high category in the survival analysis. In the epigenetic aging analysis, the alcohol index at the time of blood sampling was binned to three levels (non-drinker=0g per day, infrequent <1.3g per day, frequent ≥1.3g per day) due to the relatively small sample sizes for different categories within each reproductive class.

Education. Higher education has been shown to be associated with lower EAA, later-life morbidity and delayed mortality55,56. Education is also hypothesized to reflect differences in socio-economic status and social networks57. Lifetime years of education were compiled from answers about attained educational level in 1975 and 1981 questionnaires, with further details described in Silventoinen et al. (2004)58. We used education as a three-level (primary ≤7 years of lifetime education, secondary ≤10 years, tertiary >10 years) categorical variable in our analysis for both mortality and epigenetic aging.

Handling of missing data
Not all participants answered the 1975 questionnaire, and thus had missing values for our four covariates used as indicators for common risk factors (alcohol use, smoking, BMI and education, Supplementary Table 3). For our survival analysis sample (n=14,836), there were 10,975 women in the complete case sample, 164 women with only one missing covariate, 17 women with two missing covariates, 325 with three missing covariates and 3355 women with all four missing covariates. We evaluated the nature of missingness in R using the mcar_test from the R package naniar59, which showed our four covariates (alcohol use, smoking, BMI, and education) were not missing completely at random (MCAR; Little’s χ2(26)=94.3, p=1.13×10−9). Thus, we assumed the data were missing at random (MAR)60. For survival analysis, we used Multiple Imputation (MI) to handle missing covariate values using Mplus, assuming that these values are missing at random (alcohol use, smoking, BMI and education). MI is assumed to produce unbiased parameter estimates and standard errors, as all the available data is utilized to estimate missing variables60. We generated 20 datasets (as recommended in Little et al. 201360) with imputed covariates, and the subsequent modeling was conducted in the same way as the complete case models, by handling multiple imputed datasets with the R package mitools61. For comparison, we also ran the survival analyses with the complete case sample without imputation (Supplementary Fig. 3). For epigenetic aging data (n=1054), there were no participants with missing education information. There were 890 women in the complete case sample, 35 women with one missing covariate, 132 women with two missing covariates and 60 women with all three (alcohol use, smoking and BMI missing). We analysed epigenetic age acceleration with Bolck-Croon-Hagenaars (BCH) approach using Full Information Maximum Likelihood, which takes into account all observed data, including missing data (see section Epigenetic aging).

Epigenetic aging data
Blood sampling, DNA extraction, beadchips. Epigenetic age was assessed for a subsample of the participants (n=1054) based on peripheral blood samples taken at an age ranging from 36 to 8922. High molecular weight DNA was extracted from the blood samples with standard automated protocols and bisulfite converted using EZ-96 DNA Methylation-Gold Kit (Zymo Research, Irvine, CA, USA) according to the manufacturer’s protocol. DNA methylation levels were measured using Illumina’s Infinium HumanMethylation450 BeadChips (450k) or the Infinium MethylationEPIC BeadChips versions 1 and 2 (EPICv1 and EPICv2), which quantify methylation levels at single-nucleotide resolution of over 450,000, 850,000 and 935,000 CpG-sites, respectively. Of the samples included in the present study, 250 were assayed using the 450k array, and 667 and 200 using the EPICv1 and EPICv2 arrays. Both twins in a pair were assayed with the same platform.

Preprocessing DNA methylation data. DNA methylation data were preprocessed in R keeping data from each array type separate. Quality control and control probe-based quantile normalization were conducted using the R package meffil62. We validated our dataset by ascertaining that all methylation-based predictors of sex were indeed predicting female sex. We discarded samples where 1) median methylated signal over all CpG sites was more than 3SD from the expected, based on regression of median methylated signal by median unmethylated signal in all samples, 2) BeadChip inherent control probe values deviated over 5SD from from the overall control probe mean, 3) over 20% of probes had detection p value of over 0.01, and 4) over 20% of their probes had less than 3 detected beads. Additionally, a CpG site was removed in all samples if 1) there was only background signal in over 20% of the samples (detection p value >0.05), and 2) if the bead count of a certain probe was less than 3 in over 20% of samples. Further, all probes in sex chromosomes as well as cross-reactive probes and polymorphic CpGs were removed63,64. Raw probe intensities for each sample were adjusted to conform to its set of normalized quantiles based on 16, 15 and 16 first control probe principal components, for 450k, EPICv1 and EPICv2 platforms respectively (Supplementary Fig. 6). Beta Mixture Quantile Normalization65 to adjust the beta-values of type II design probes into a statistical distribution characteristic of type I probes was conducted using the R package wateRmelon66. Methylation beta values were obtained by dividing the intensity of methylated sites (M) by the sum of intensities from methylated and unmethylated sites (Beta=M/(M+U) where U represents the intensity of unmethylated probes). These processed beta values were then used as input for calculating epigenetic aging.

Epigenetic aging data. Biological age was determined from DNA methylation beta values using three published algorithms: principal component versions17 of GrimAge16 and PhenoAge15 (referred to as PCGrimAge and PCPhenoAge), as well as DunedinPACE14. We used only the probes that were present in all platforms after preprocessing (for PC-clocks 66,260 CpGs/78,464 CpGs; for DunedinPACE 144 CpGs/173 CpGs). Epigenetic ages of the above algorithms were calculated in R using the scripts provided in the corresponding publications. Epigenetic age acceleration was determined as the residual of epigenetic age linearly regressed against an interaction between chronological age and platform, and bisulfite-conversion plate in our sample of 1,054 participants.

Mortality data
The mortality follow-up ended on 31 December 2020, when dates of death were retrieved from the Digital and Population Data Services Agency Finland. The survival time for each woman was defined in years until the date of emigration, death, or end of follow-up, whichever came first.

Statistical analysis
Identifying reproductive trajectories by latent class analysis. Reproductive trajectories were identified from lifelong live birth data of 10,783 women using latent class analysis (LCA)9 with the Mplus software, analysis launched via R using the package MplusAutomation66. The process iterates from random starting values to explore different regions of the parameter space and find values that maximize the likelihood of observing the data given the model. The estimated parameters include class-specific item response probabilities (to give birth) for each indicator (ages at childbirth, see below), as well as class probability parameters, which are used to calculate individual-level posterior probabilities of belonging to each class. The estimation was done by using maximum likelihood estimators with robust standard errors (MLR) using a sandwich estimator (MLR type=COMPLEX)67. The sandwich estimator was used as it is robust to violations of non-independence among observations resulting from shared inheritance. For initiating the estimation we used 500 random sets of parameter values and 20 final stage optimization rounds; the number of starting values and optimizations was increased when needed for the model to converge to global maximum68. Global maximum was assumed when the final models converged to the same likelihood at the final stage optimization rounds68.

The optimal model solution was determined following guidelines by Sinha et al. (2021)69, suggesting to select the model with the smallest number of classes that best fits the data. Fit was assessed with Akaike’s information criterion (AIC)70 and the sample size adjusted Bayesian information criterion (saBIC)71,72, which both assess model accuracy based on likelihood but penalize overfitting. If these criteria disagreed we prioritized saBIC as it is recommended especially for analyses where sample size is large73,74. After prioritizing information criterion, we aimed to maximize the relative sample size of the smallest class, entropy (a measure of class separation) and class-specific posterior probabilities. We tested models with the number of live births summed over 1–4 year bins as the indicator variables (for example in the 3-year age bins, all births younger than 18 years, 18–20 years, 21–23 years,…, older than 44 years). For each age bin, we fitted latent class models with 1 to 8 classes. The number of live births within 3-year age bins were selected as indicator variables for the final latent class analysis, with 6 latent classes. Comparisons between all the tested models can be found in Supplementary Table 1. In the final model solution, we added women with no reported live births i.e. the nulliparous as an additional class, for which an individual was assigned a probability of 0.9999999 if they had no live births (and 0.0000001 as the class-specific posterior probabilities), and 0.0000001, if they had given birth during the study period and were included in the LCA solution.

Parametric bootstrapping of the latent class solution. The robustness of the final LCA model was tested with parametric bootstrapping in R. We simulated 100 datasets consisting of 10,783 individuals. For each age bin i, the individual probability of having n live births was calculated as the product of the posterior class belonging probability matrix (I classes × S individuals), by the latent class matrix (N possible live births × I classes), yielding a class-specific (N possible live births × S individuals) probability matrix. For each bootstrap replicate, an actual number of live births in each age class for each individual was drawn randomly according to these probability matrices. For each simulated dataset, the latent class model with 6 classes was then fitted.

Assessing differences in survival between the reproductive classes. The differences in survival between the latent classes including the nulliparous class were modeled in R. To account for the non-independence between women from the same twin pair, we conducted a mixed Cox proportional hazards model, where twin pair ID was used as a random effect using R packages survival75 and coxme76. To reflect latent class admixture proportions, we used proportional assignment: each individual was assigned to all of the 7 classes (6 LCA classes and the nulliparous) simultaneously with weights equal to class-specific posterior probabilities77. We used class 3 as the reference class in our survival analyses to reflect a reproductive pattern that is common within this population, as this class had the highest probability of childbirth at age 27–29, while the average age at childbirth was 27.3 in our total sample (Supplementary Results).

The follow-up started at birth, and as the study participants only entered the study provided they answered the 1975 questionnaire, we used 1975 as an entry date for left truncation in our model. Right-censored date of exit was either death, emigration or end of follow-up, whichever came first. Two models were fitted: an unadjusted model included only historical birth cohort and twin pair ID variables, whereas an adjusted model additionally included smoking, alcohol use, and BMI in 1975, and lifetime education (see variable description above). Statistically significant differences in mortality were defined by hazard ratios along with their 95% confidence intervals from the Cox proportional hazards model.

Assessing differences in epigenetic aging between the reproductive classes. The differences in epigenetic age acceleration between the 6 latent reproductive classes and the nulliparous women were modeled in MPlus using the Bolck-Croon-Hagenaars (BCH) approach which controls for measurement error in the classification78. We used class-specific weights as training data to model the association between epigenetic age acceleration and latent classes. Similarly to survival analysis, we constructed two different models: an unadjusted model, which included independent variables for the twin pair ID and chronological age, and the adjusted model, which additionally included variables for smoking, alcohol use, BMI at the time of blood sampling, and lifetime education.

Assessing effects of parity and age at first child on later health. The effects of the number of live births and age at first child were explored using categorical child count with five levels (one, two, three, four and five or more birth givings), and age at first child, respectively. All the models were run without imputation of missing covariates, and thus used the complete case sample. First, we modeled survival with mixed Cox proportional hazards models with adjustments for left truncation, alcohol use, smoking, BMI and education, using the same modeling methods in R as for the main analysis with reproductive classes. Women with three live births were used as a reference group for the survival analysis. Second, we modeled epigenetic age acceleration with linear mixed effects models in R using the package lme479. These models were adjusted with chronological age, alcohol use and tobacco use, and BMI at the time of blood sampling and lifetime education. All of the above models included twin pair ID as a random effect to account for related individuals. The epigenetic aging estimates are reported using estimated marginal means with the R package emmeans80.

Alternative analysis strategies. To test whether our results would be impacted by the modeling framework, we modeled our original Cox proportional hazards models in R with cluster-robust standard errors (Supplementary Table 5). The models were otherwise the same as our original models, but twin pair ID was used as a clustering variable instead of a random intercept. Finally, we down-sampled our dataset to a single individual per family, removing the twin structure entirely from our data. This singleton approach was repeated 100 times, each time randomly selecting one individual per family. For each random singleton subset, we modeled our left-truncation adjusted Cox proportional hazard models weighted with class-specific posterior probabilities in R, as well as our epigenetic age acceleration (EAA) models using the Bolck-Croon-Hagenaars method in Mplus, following the same adjustment methods as in our main analyses (excluding random intercept). We summarized the hazard ratios, estimates and 95% confidence intervals across these models in a fixed-effects meta-analysis in R using the package metafor81.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The revision has addressed the concerns raised by the referees to the extent that the structure of the data set allows, and the limitations of the data set and analysis are now more clearly communicated.", "reviewer_2": "The authors have addressed my comment. I note that the approach the authors opt for is vulnerable to bias in the case that unmeasured factors are correlated with both the outcome and the unit of the random effect (in this case, the family unit). Such factors are almost certainly present -- genetics, unmeasured variability in socioeconomic circumstances, family norms in health behaviors, etc. There is a formal test to determine if random effects are valid in this setting proposed by Houseman (e.g. see here https://www.sciencedirect.com/topics/social-sciences/hausman-test). It almost always rejects the random effects approach, though. An alternative might be to evaluate the consistency of results using alternative methods to account for the non-independence of observations within families: clustering standard errors at the family level, generalized estimating equation (GEE) models. At minimum, the authors should acknowledge the vulnerability of their analysis to this type of confounding.", "reviewer_3": "The authors have voluminously responded to previous reviewer comments. It is strongly suggested to just put these limitations into the manuscript rather than inching up to them, writing 4 paragraphs in response, and then disregarding in the text.\nSome points for perspective:\nPregnancy is a well-known physiological stressor. As noted by multiple reviewers, this isn’t a new finding or surprising. Well documented phenomena across populations.\nThe sub-sample of epigenetic measures taken between 1994-2020, along with the table 1 data on the underlying population could be a case-study for selection bias, and plausible avenues for why what is observed in the data is observed. Based upon the voluminous responses to the reviewers, consulting a scientist who understands this concept would have been a better use of everyone’s time. In short, there appears to be survival based requirement to be included, along with a response requirement. The reader doesn’t know the baseline denominator of potential subjects for the baseline assessment. The epigenetic data was a convenience sample and not representative of the overall cohort.\nWhile the goal of the research may indeed have a causal lens, the methods applied are not estimating a causal effect, but the language in the text throughout implies causality. This was not changed, despite reviewer comments.\nThe Cox regression survival analysis does not define how follow up was calculated or what time 0 was, making it impossible to interpret. It takes a puzzle master to figure out when the various covariates were measured/included as well, and still unclear how they align.\nBMI is not a lifestyle factor.\nYes, the application of latent class modeling to derive “distinct” classes helps statistically, but what any of the classes represent beyond the nulliparous and those with a high number of pregnancies. Despite this, there’s nothing really going on except on the tails of the distribution, and despite nulliparous women being at the highest risk compared to whatever latent class 3 is, the text largely ignores this and spends most of the time speculating and interweaving supported theories with evolutionary reproduction factors and focusing the conclusions on this.", "reviewer_4": "The authors have been very responsive to the reviews and revised the manuscript satisfactorily."}, "author_response": "We thank the reviewer for bringing this issue to our attention - indeed this needed to be clarified better in our text. We agree with the reviewer that random-effect models indeed impose a “restricted” distribution of shared family-level intercepts, and that this can under-estimate or in extreme cases bias estimates.\nA within-pair-difference model, as suggested here, would indeed be able to do that - although it would raise other questions about covariate formulation, as the modelling unit becomes the twin pair and not the individual, and consequently we need to account for both average value and pair-difference in each included covariate (the latter having no standard definition for e.g. study level, quantitative tobacco use, etc). And most importantly, it would strongly reduce our sample size, as we only focus on women here, and we would not be able to include the female co-twin of opposite-sex dizygotic pairs in our model anymore.\nOur approach using family level random effect is thus a compromise that aims at correcting for family relatedness, where the introduction of the random effect is only meant to avoid violating the assumption of independence of observations that we would encounter when using a standard linear regression. Our approach has been described with regards to twin models by Fasola et al. (2021) as generalized linear mixed models (GLMMs), where the authors have also compared the approach to other methods used in twin analysis, such as the within-pair model. GLMMs are commonly used in twin data analyses, as they allow a response variable to be explained by both fixed effects (which apply uniformly to all individuals) and random effects (which vary by family, capturing sources of similarity between twins), enhancing the robustness of statistical inferences by appropriately modeling the dependencies inherent in twin data. This method has been widely used in twin studies (Yuan et al. 2014, Fasola et al. 2021), also in our cohort (Huang et al. 2020, Gustavson et al. 2025).\nHowever, we understand that describing the analyses as \"twin model analyses\" in the manuscript may be misleading and we agree that there was a lack of clarity about how the models and their results should be interpreted, which we have now corrected in the manuscript. We now clearly state that the family ID is only meant to correct for statistical non-independence between the twins from the same family throughout the manuscript. In addition, we have made sure that there is no confusion regarding potential corrections for shared factors between co-twins, as our models do not allow this.\nWe agree that pregnancy, as well as the associated burden of raising children are well-known physiological stress factors, which was the basis for initiating this study.\nIn the Finnish Twin Cohorts the high response rates in follow-ups since the founding of the study in 1975 indicate the FTC is broadly representative of the general population (Kaprio et al. 2019). Long-term analyses show no significant differences from the general population in mortality, cancer incidence, or genetic studies (Kaprio 2013, Skytthe et al, 2019). For our study, we utilized all available DNA methylation data among the cohort participants collected from various targeted studies and general expansions without phenotype-based selection. Studies performed utilising epigenetic clocks in the cohort data show very consistent results with other population studies. However, we acknowledge that there might still be some selection bias related to the prerequisite to be alive and answer the questionnaire in 1975 and that some epigenetic data has been gathered in purpose to target specific traits, and we have now added this to the Limitations-section of our Main text:\nThird, while the Finnish Twin Cohort has high response rates and has been shown to be broadly representative of the general population7,42,43, we acknowledge that some selection bias may be present due to the requirement for participants to be alive and respond to the baseline survey in 1975. Additionally, the available DNA methylation data was collected from various targeted studies and general expansions, meaning that some samples were selected to investigate specific traits rather than being fully representative of the overall cohort. Additionally, the available DNA methylation data were drawn from a combination of targeted studies and broader cohort expansions. As a result, some samples were specifically selected to investigate particular traits, rather than to represent the full cohort. Consequently, analyses based on this subsample of epigenetic data should be interpreted with caution, as the participants included may be subject to selection bias and may not fully reflect the general population by their physiological condition.\nWe have now clarified our manuscript to indicate that our findings are associations, and that any discussion about causality is speculative. We have also added this point to our limitations-section:\nFinally, while our study examines associations between reproductive history, epigenetic aging, and survival, the observational nature of our data and the methods used do not allow for causal inference. Although we employ rigorous statistical approaches to control for potential confounders, our findings should be interpreted as associations rather than causal effects. Future research using for example, quasi-experimental designs, such as Mendelian randomization or other instrumental variable approaches, may help to better establish causal relationships.\nWe have now specified the start of the follow-up and the adjustment for delayed entry throughout the manuscript, and described our method in the Main text and the Supplementary methods (see below). We have also specified the information on when the covariates were collected for the survival and epigenetic aging analyses.\nIn the Main text: Methods line 537:\nThe start of the follow-up was birth, and we used left truncation to account for delayed entry from the time of the first questionnaire in 1975.\nSupplementary methods: 3.2 Survival:\nThe follow-up started at birth, and as the study participants only entered the study provided they answered the 1975 questionnaire, we used 1975 as an entry date for left truncation in our model.\nWe have now changed the description of our “lifestyle” adjustments to “risk factors” throughout the manuscript.\nWe appreciate the reviewer’s comment and agree that the most pronounced differences in survival and aging outcomes appear at the extremes of the reproductive distribution, particularly for nulliparous women and those with a high number of pregnancies. However, our latent class modeling approach was designed to explore whether there were additional meaningful distinctions in survival and aging patterns beyond these well-known extremes. Since we could not assume in advance what the result would be, our goal was to identify potential heterogeneity within the typical reproductive patterns (classes 2–5) rather than focusing solely on predefined groups.\nOur findings indicate that, after adjusting for confounders, differences within these middle-range classes were minimal. While this outcome may suggest that reproductive timing and moderate variations in offspring number are linked to smaller variations with our outcomes than initially expected, our approach was necessary to establish this empirically. To further support this conclusion, we also conducted analyses using the number of births and age at first child as continuous predictors (see Extended Data Figure 7 and Supplementary Results), reinforcing our confidence that no major overlooked patterns exist, particularly concerning childbirth timing. We have described this rationale in our Main text under the subheading: Typical reproduction is the optimal reproduction."}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have addressed my comment. I note that the approach the authors opt for is vulnerable to bias in the case that unmeasured factors are correlated with both the outcome and the unit of the random effect (in this case, the family unit). Such factors are almost certainly present -- genetics, unmeasured variability in socioeconomic circumstances, family norms in health behaviors, etc. There is a formal test to determine if random effects are valid in this setting proposed by Houseman (e.g. see here https://www.sciencedirect.com/topics/social-sciences/hausman-test). It almost always rejects the random effects approach, though. An alternative might be to evaluate the consistency of results using alternative methods to account for the non-independence of observations within families: clustering standard errors at the family level, generalized estimating equation (GEE) models. At minimum, the authors should acknowledge the vulnerability of their analysis to this type of confounding.", "reviewer_3": "The authors have largely responded to everything and updated the manuscript. One point remains which needs updating, R 3.2:\nThe following statement should be removed from the text: \"Nevertheless, studies utilizing epigenetic clocks within this dataset have produced results that align well with other population studies44-46, suggesting that major selection effects are unlikely.\"\nThis is incorrect and not an approach or valid way to address or assuage the selection bias."}, "author_response": "We thank the reviewer for raising the important issue of unobserved heterogeneity—specifically, the assumption in mixed-effects models that the random effect (in our case family-level intercept) is not correlated with the predictors (such as reproductive class or covariates). We agree that if this assumption is violated, it could bias the model estimates and their interpretation.\nIn response to this concern, we conducted several additional analyses to assess the robustness of our model and the potential for bias. These results are developed below and have been added to the Main text and Supplementary Methods.\nDue to the large size of our dataset (~9,500 families), the fixed-effect model necessary for the Hausman test would not converge on the Finnish IT center for science computing grid. As a result, we proceeded directly to testing generalized estimating equation (GEE) models as suggested by the reviewer. As we are not aware of any GEE implementation for Cox proportional hazards models matching our data structure (left-truncated, right-censored data, proportional weighting by latent class and time-to-event outcomes), we employed Accelerated Failure Time (AFT) models using the GEE framework, which allow modeling of survival time (rather than hazard ratios) while adjusting for right censoring, latent class weights, and covariates. However, the only currently available implementation in R to our knowledge (aftgee) is unable to account for left truncation (i.e. the fact that participants entered the study in 1975 at different ages). The results of the AFT-GEE model adjusted only for historical cohorts supports our original conclusions (see results below). In this “unadjusted” model, the nulliparous women have a time ratio of 0.95 (95% CI: 0.94–0.96), indicating a 5% shorter survival time, while class 1 and 2 have time ratios of 0.98 (CI: 0.97–0.99), corresponding to 2% reduction in survival time. Finally, class 6 has a time ratio of 0.97 (CI: 0.96–0.99), corresponding to a 3% reduction in lifespan. Unfortunately, as in the case of the Hausman test, the risk factor (alcohol use, smoking, BMI and education) adjusted model did not converge in reasonable time (10 days) on our computing grid, despite repeated attempts.\nBecause of this difficulty, and of AFT-GEE models’ inability to account for left truncation, we report this modeling strategy here for completeness but do not include it in the main manuscript.\nAs suggested by the reviewer, we also re-estimated our survival models without random effects but with cluster (here family) robust standard errors (Extended Data Table 5, reproduced below). The survival model estimates with cluster robust standard errors are consistent with our original results, providing additional evidence for the robustness of our findings. For epigenetic age acceleration analyses, the models presented in the manuscript already use cluster robust standard errors Bolck–Croon–Hagenaars method with type=COMPLEX in Mplus).\nHere, we propose a third and more straightforward method to test for the effect of unmeasured family-level effects on our results: as a control, we simply down-sampled our dataset to a single individual per family, removing the twin structure entirely from our data. This “singleton” approach was repeated 100 times, each time randomly selecting one individual per family. For each random singleton subset, we modeled our left-truncation adjusted Cox proportional hazard models weighted with class-specific posterior probabilities in R, as well as our epigenetic age acceleration (EAA) models using the Bock-Croon-Hagenaars method in Mplus. We summarized the hazard ratios, estimates and 95% confidence intervals across these models using the function rma(method=\"FE\") from the R package metafor(). Comparisons of hazard ratios (survival) and estimates (EAA) along with their confidence intervals from these singleton models to our original models are presented in Extended Data Table 6, and the estimates from the individual models are shown in Extended Data Figure 8 (reproduced below). Sample sizes vary slightly in non-imputed Model 2 due to differences in the number of individuals with missing data across the random subsets. In contrast, sample sizes remain consistent for Survival Model 1 (complete-case model which includes no missing data) and for the imputed Model 2 (where missing values are replaced by imputed values). Sample sizes remain consistent in the non-imputed BCH Model 2 because it employs Full Information Maximum Likelihood (FIML), which retains all observations by using available data without excluding cases due to missing covariate values. These singleton analyses unambiguously support our original conclusions regarding the between-class differences in both survival and EAA.\nWe have now incorporated these results to the manuscript as follows:\nMain text:\nLine 181: We used alternative modeling strategies to assess the robustness of our results. First, we reanalyzed survival time using fixed-effects models with cluster-robust standard errors. Second, we implemented a “singleton” approach by drawing 100 random subsamples, each including one individual per family from the 9,496 families, and reran our survival models on these samples. Both approaches consistently showed elevated hazard ratios for nulliparous women, early mothers, and the high lifetime reproductive group (classes 0, 1 and 6). The increased risks for nulliparous women and the high lifetime reproductive group remained statistically significant after adjusting for risk factors, underscoring the robustness of these associations (Supplementary Methods, Extended Data Fig. 8, Extended Data Tables 5 and 6).\nLine 221: These results were also supported by our “singleton” approach with 100 random draws of individuals from each of the 611 families in our epigenetic age acceleration analyses (Extended Data Table 6, Extended Data Fig. 8).\nSupplementary methods:\n3.5 Alternative modeling strategies\nTo test whether our results would be impacted by the modeling framework, we modeled our original Cox proportional hazards models with cluster robust standard errors (Extended Data Table 5). The models were otherwise the same as our original models, but Family ID was used as a clustering variable instead of a random intercept. Finally, we down-sampled our dataset to a single individual per family, removing the twin structure entirely from our data. This “singleton” approach was repeated 100 times, each time randomly selecting one individual per family. For each random singleton subset, we modeled our left-truncation adjusted Cox proportional hazard models weighted with class-specific posterior probabilities in R, as well as our epigenetic age acceleration (EAA) models using the Bock-Croon-Hagenaars method in Mplus. We summarized the hazard ratios, estimates and 95% confidence intervals across these models in a fixed-effects meta-analysis using the R package metafor44.\nWe thank the reviewer for noting this. We have now removed the sentence mentioned by the reviewer, and modified the limitations section at line 333 onwards:\n“Additionally, the available DNA methylation data were drawn from a combination of targeted studies and broader cohort expansions. As a result, some samples were specifically selected to investigate particular traits, rather than to represent the full cohort. Consequently, analyses based on this subsample of epigenetic data should be interpreted with caution, as the participants included may be subject to selection bias and may not fully reflect the general population by their physiological condition.”"}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors have addressed my comments", "reviewer_3": "The authors have sufficiently addressed the review comments this reviewer raised."}, "author_response": ""}]

Decision letter