Large-scale whole-genome sequencing reveals the landscape and health implications of de novo mutations.

DNMs are spontaneous genetic alterations that arise within the current generation and are absent from parental genomes. Although neutral or beneficial DNMs contribute to evolution and genetic diversity, those disrupting functional genomic elements may result in severe phenotypic consequences, such as rare sporadic malformation syndromes 1 , 2 . Over the past decades, considerable efforts have been made to assess the rate and etiology of parent-of-origin DNMs in humans, which account for nearly 90% of all DNMs 3 , 4 . Parental age at conception is a well-established determinant of germline DNMs, with nearly three-quarters of age-related DNM accumulation arising from the paternal germline 4 – 9 . Epidemiological studies have also consistently associated advanced parental age with adverse birth and developmental outcomes 10 , 11 , leading to the hypothesis that age-related DNM accumulation contributes to these effects 12 . However, this hypothesis has not been directly evaluated within a single prospective cohort integrating genomic and longitudinal follow-up data. Concurrently, the rise in delayed parenthood has been accompanied by the increasing use of ART. Our previous study suggested that ART may be associated with altered DNM burden in offspring 13 . However, ART encompasses multiple clinical procedures, making it unclear whether the observed DNM changes arise from parental characteristics, sperm retrieval and ovarian stimulation procedures, embryo culture conditions or the bypassing of natural selection. Moreover, whether ART-associated mutational changes translate into measurable effects on offspring health remains unclear. Beyond parent-of-origin DNMs, post-zygotic mutations that occur during early embryonic development 14 – 16 can also influence offspring phenotypes 17 , yet their relationship with ART has not been systematically investigated. In the present study, we performed whole-genome sequencing (WGS) of 24,030 individuals from 7,851 parent−offspring families in the Jiangsu Birth Cohort (JBC). We systematically investigated how parental age and ART influence the rate and spectrum of parent-of-origin DNMs and early post-zygotic mosaic mutations (EPZMs), whether specific ART procedures exert distinct mutational effects and whether these mutations statistically account for associations among parental age, ART use and offspring birth or developmental outcomes (Fig. 1a ). Fig. 1 Study workflow and the mutational spectrum of DNMs. a , Graphical abstract of the study: WGS of 8,328 infants from 7,851 parent−offspring families identified approximately 0.42 million high-quality DNMs, of which 10.22% were EPZMs and 89.78% were parent-of-origin DNMs. Socio-demographic information and reproductive histories of the recruited couples from the JBC study, along with detailed ART-related characteristics, were collected using standardized questionnaires and medical records (left). We characterized the mutational spectra of both parental DNMs and EPZMs (middle) and evaluated their associations with birth outcomes, congenital anomalies and neurodevelopmental performance at 1 year of age (right). Proportions of paternal and maternal DNMs were calculated by first excluding EPZMs and then presenting the ratio of paternal and maternal origin among the remaining mutations based on the inferable parent-of-origin assignments. Selected icons illustrating ART procedures in this figure were adapted from Servier Medical Art ( https://smart.servier.com ), licensed under CC BY 4.0 ( https://creativecommons.org/licenses/by/4.0/ ). b , Relative frequencies of mutational classes of parent-of-origin and EPZM dnSNVs. c , Relative frequencies of trinucleotide mutational contexts in parent-of-origin and EPZM dnSNVs (top). Heatmap of mutation profile differences split by nucleotide substitutions and context (bottom). Blue indicates overrepresentation in parent-of-origin or maternal dnSNVs, respectively, and red indicates overrepresentation in EPZM or paternal dnSNVs, respectively. d , Relative frequencies of parent-of-origin and EPZM dnSNVs across different genomic regions. For b − d , proportions were compared using two-sided Fisherʼs exact tests. Bonferroni correction was applied. * P  < 0.05; ** P  < 0.01; *** P  < 0.001. Exact P values are provided in the Source Data file. ~, approximately; NS, not significant; OR, odds ratio. Source data a , Graphical abstract of the study: WGS of 8,328 infants from 7,851 parent−offspring families identified approximately 0.42 million high-quality DNMs, of which 10.22% were EPZMs and 89.78% were parent-of-origin DNMs. Socio-demographic information and reproductive histories of the recruited couples from the JBC study, along with detailed ART-related characteristics, were collected using standardized questionnaires and medical records (left). We characterized the mutational spectra of both parental DNMs and EPZMs (middle) and evaluated their associations with birth outcomes, congenital anomalies and neurodevelopmental performance at 1 year of age (right). Proportions of paternal and maternal DNMs were calculated by first excluding EPZMs and then presenting the ratio of paternal and maternal origin among the remaining mutations based on the inferable parent-of-origin assignments. Selected icons illustrating ART procedures in this figure were adapted from Servier Medical Art ( https://smart.servier.com ), licensed under CC BY 4.0 ( https://creativecommons.org/licenses/by/4.0/ ). b , Relative frequencies of mutational classes of parent-of-origin and EPZM dnSNVs. c , Relative frequencies of trinucleotide mutational contexts in parent-of-origin and EPZM dnSNVs (top). Heatmap of mutation profile differences split by nucleotide substitutions and context (bottom). Blue indicates overrepresentation in parent-of-origin or maternal dnSNVs, respectively, and red indicates overrepresentation in EPZM or paternal dnSNVs, respectively. d , Relative frequencies of parent-of-origin and EPZM dnSNVs across different genomic regions. For b − d , proportions were compared using two-sided Fisherʼs exact tests. Bonferroni correction was applied. * P  < 0.05; ** P  < 0.01; *** P  < 0.001. Exact P values are provided in the Source Data file. ~, approximately; NS, not significant; OR, odds ratio. Source data

The present study was conducted within the JBC, a family-based prospective cohort with longitudinal follow-up and ongoing enrollment of families receiving ART treatment and families with spontaneous pregnancies at clinics in Jiangsu, China. Couples undergoing ART were enrolled at obstetrics clinics before treatment, whereas couples with spontaneous conceptions were recruited in the first trimester (8−14 weeks of gestation) when establishing health records at clinics. All ART and natural conceptions were followed with the same protocol throughout pregnancy and postnatally. The overall design and overarching goal of JBC are summarized here and are described in detail elsewhere 51 . In the present study, participants were recruited from the Women’s Hospital Affiliated to Nanjing Medical University, Suzhou Hospital Affiliated to Nanjing Medical University, Changzhou Maternity and Child Health Care Hospital Affiliated to Nanjing Medical University and Wuxi Maternity and Child Health Care Hospital. Blood samples were collected from both biological parents at recruitment, and fingertip blood samples were obtained from the offspring at the 1-year follow-up visit. Socio-demographic information (for example, age, sex, residence area, education and annual household income), pre-pregnancy weight and height and reproductive history of the recruited couples were collected with standardized and structured questionnaires by face-to-face interviews. Clinical information, including detailed ART-related characteristics, pregnancy complications (for example, gestational diabetes and hypertensive disorders during pregnancy) and neonatal outcomes, was extracted from medical records of the participating hospitals. The JBC biospecimen banking and genomic study was conducted in accordance with relevant ethical regulations and approved by the ethics committees of Nanjing Medical University and the participating hospitals. All participants provided written informed consent and received no financial compensation. Inclusion criteria were as follows: (1) couples delivered live births; (2) both spouses completed the standardized and structured questionnaires; (3) infants completed the 1-year follow-up and health assessment; and (4) both parents donated blood samples at recruitment, and their infants donated a fingertip blood sample at the 1-year follow-up. A total of 24,318 participants from 7,947 parent−offspring trios met the above criteria and were subjected to WGS (Fig. 2 ). ART was defined as conception through IVF with or without ICSI. Parental infertility diagnoses, specific ART procedures and medication usage were also exposures of interest. Infertility diagnoses of ART-conceived couples included polycystic ovary syndrome, ovarian dysfunction, endometriosis, pelvic or tubal factors and OAT. Couples with more than one infertility diagnosis were included in each subgroup analysis. ART procedures included embryo transfer cycle type (fresh/frozen), stage of embryo transfer (cleavage/blastocyst), fertilization method (IVF/ICSI) and controlled ovarian stimulation protocols. Medication usage during ART procedures included trigger strategies, subtypes of follicle-stimulating hormone (FSH) supplements and drug supplements of GnRH-ant, luteinizing hormone and clomiphene. We were authorized to extract health outcome data for cohort participants, including infant sex, gestational duration, birth weight and birth defects, from the electronic medical records (EMRs) of the participating hospitals where women underwent antenatal examinations and gave birth. Preterm birth was defined as any live birth occurring before 37 weeks of gestation. Low birth weight was defined as birth weight less than 2,500 g. Birth defects were diagnosed by licensed clinicians, and the diagnoses were coded using the International Classification of Diseases, 10th edition (ICD-10). Detailed diagnoses and corresponding ICD-10 codes during prenatal follow-ups, at birth and in early infancy were obtained from the EMR. Any new diagnoses of birth defects were identified through telephone follow-ups with parents at 42 days and 6 months postpartum. All cohort children received health examinations at the child healthcare department at 1 year of age. Children with parent-reported birth defects or those suspected of having defects after a systematic physical examination were referred to licensed specialists for definitive diagnoses. Birth defects were categorized into 11 subtypes according to ICD-10 codes, including nervous system defects; eye, ear, face and neck defects; circulatory system defects; respiratory system defects; cleft lip and cleft palate defects; digestive system defects; genital organ defects; urinary system defects; musculoskeletal system defects; chromosomal abnormalities; and other birth defects 52 . The Bayley Scales of Infant and Toddler Development, 3rd edition, screening test (Bayley-III Screening Test) was used to evaluate the neurodevelopment of infants aged 11−12.5 months in this study, which consists of five domains: cognition, receptive communication, expressive communication, fine motor and gross motor. Using the Illumina HiSeq or NovaSeq platforms, 150-bp paired-end WGS was performed according to the manufacturerʼs instructions (Illumina), with a target depth of approximately 30× (Supplementary Fig. 11 ). Of the 7,851 sequenced parent−offspring families, 331 families (4.22%) were sequenced on the Illumina HiSeq platform, and the remaining 7,520 families (95.78%) were sequenced on the NovaSeq platform. Raw paired-end WGS reads were subjected to adapter and low-quality sequence trimming using Trimmomatic (version 0.36) with the default parameters 53 . The FastQC package ( http://www.bioinformatics.babraham.ac.uk/projects/fastqc ) was used to assess the quality score distribution of the sequencing reads. Read sequences were mapped to the human reference genome (GRCh37) using the Burrows−Wheeler Aligner (BWA-MEM version 0.7.15) 54 with default parameters, and duplicates were marked and discarded using SAMBLASTER (version 0.1.24) 55 . The quality metrics (including coverage, median insert size and percentage of chimeric reads) and contamination estimation of ‘BAM’ files were obtained by Picard (version 1.70) ( http://broadinstitute.github.io/picard ) and VerifyBamID2 (version 1.04) ( https://github.com/Griffan/VerifyBamID ) 56 , respectively. The local realignment and recalibration of the reads aligned by BWA were conducted using the Genome Analysis Toolkit (GATK version 3.8) 57 . Haplotype calling was performed using HaplotypeCaller in GATK version 3.8 in GVCF mode according to GATK Best Practices. Individual gVCF files were jointly genotyped for high-confidence alleles using GenotypeGVCFs for all autosomes and the X chromosome. To estimate the genetic relatedness and sex of the study populations, we generated a temporary list of high-confidence biallelic variants with (1) a call rate greater than 99% and (2) an allele frequency greater than 1% in the study cohort. Pairwise relatedness was estimated by KING 58 , and the chromosomal sex of samples was inferred based on the inbreeding coefficient ( F ) for common variants on chromosome X with pseudo-autosomal regions excluded. Population structure was estimated by using principal component analysis in EIGENSOFT (version 2.7.1) with linkage disequilibrium−pruned autosomal biallelic variant sets ( r 2  = 0.1). Individuals were excluded if they met any of the following criteria: (1) average coverage <10×; (2) median insert size 0.05); (4) high contamination: freemix >0.05; (5) unmatched relatedness; or (6) inconsistent chromosomal sex and questionnaire information. A rigorous quality control process involving the combination of hard filters and a random forest model was implemented for variant filtering. First, variants meeting any of the following criteria were excluded: (1) with an excess of heterozygotes, indicated by an inbreeding coefficient of less than −0.3; (2) spanning deletion; (3) with variant lengths exceeding 50 bp; or (4) with call rates below 95%. Then, a random forest model was applied to distinguish true variants from potential artifacts. The random forest model was optimized and trained with the Python auto-sklearn package 59 to distinguish true variants from potential sequencing artifacts. The model operated at the allele level and jointly considers SNVs and indels. A total of 61 features were incorporated (Supplementary Table 17 ), including (1) site-level annotations generated by GATK HaplotypeCaller (for example, genotype quality, inbreeding coefficient and ReadPosRankSum); (2) allele-level annotations derived from in-house scripts (for example, variant allele fraction and allele-level quality by depth); and (3) regional annotations reflecting local genomic context (for example, GC content and presence within simple repeat sequences). Positive and negative training datasets were constructed from WGS data. The positive training set included alleles previously genotyped or confidently identified in three widely used reference resources: Omni 2.5, Mills and the 1000 Genomes Project. Because variants in these databases are predominantly common, we additionally included singleton variants observed in unrelated samples that were transmitted in at least one trio, thereby improving the modelʼs ability to distinguish rare true variants from technical artifacts. The negative training set consisted of alleles that failed standard GATK hard filters (QD  60 or MQ < 30). To further evaluate model performance, chromosome 20 was held out entirely during training, and performance was assessed on this withheld chromosome, yielding an area under the receiver operating characteristic curve (AUC) greater than 0.99. Because callable regions varied across families, we defined the consensus callable regions to minimize bias when comparing the DNM rate. We first excluded difficult-to-map regions defined by the Genome in a Bottle (GIAB) benchmark set for polymerse chain reaction (PCR)-amplified short-read sequencing data 60 , including tandem repeats, perfect or imperfect homopolymers >10 bp, segmental duplications and regions homologous to the hs37d5 decoy or ALT loci. We further excluded segmental duplication and simple tandem repeats annotated from the UCSC genome browser ( http://genome.ucsc.edu/ ). After these filters, we obtained a set of high-confidence mappable regions (approximately 2.55 Gb). Within these regions, callable bases were defined as genomic sequences with genotype quality ≥20 and a minimum sequencing depth ≥10 in each individual. The callable region for each trio was calculated as the intersection of callable bases across all three individuals. On average, the callable region accounted for approximately 87.23% (interquartile range (IQR): 86.32−87.95%) of the high-confidence mappable genome. Autosomal SNVs that passed the quality control process were subjected to DNM analysis with DeNovoGear (version 1.1.1) for each parent−child trio. The identified DNMs were further filtered with the following criteria: (1) minor allele frequency = 0 in all parents included in this study cohort; (2) the offspring must be an alternative allele carrier; (3) with a minimum coverage ≥10 for parents and offspring; (4) no more than one read supporting the alternative allele in the child’s parents; (5) maximum parental allelic fraction of 0.05; (6) minimum offspring allelic fraction of 0.20; and (7) were located in the callable regions. To further eliminate potential mapping artifacts and paralogous alignments, we extracted the flanking sequences of candidate variants and performed local realignment using BLAT. Variants were removed if they exhibited evidence of multi-mapping to paralogous regions or poor alignment quality, defined as BLAT GScore <50, alignment score 10% or Tgap >50. Variants showing significant strand bias were also excluded. To account for the different read depths across infants, the standardized number of DNMs was estimated by dividing the variant number by the callable base number and then scaling to the whole genome size of 2.7 × 10 9  bp. To compare the relative frequencies of specific mutation types of DNMs across different groups, the DNMs set was grouped into six classes according to nucleotide substitution (that is, C > A, C > G, C > T, T > A, T > C and T > G). C > T substitutions were further divided into two classes based on the CpG context (that is, C > T on CpG sites and non-CpG sites). Of these DNMs, variants with VAF less than 0.35 were defined as EPZMs 61 , 62 . We evaluated the allelic fraction threshold of 0.35 for defining EPZMs in a publicly available database (CEPH/Utah pedigrees study), which included DNMs validated by three-generation pedigrees 15 . The threshold of an allelic fraction of 0.35 captured 75% (303/404) of post-zygotic mutations while limiting the false-positive rate for germline variants to 7.54% (287/3,807) (Supplementary Fig. 12 ), consistent with the previous findings 16 . We also assessed the concordance of EPZMs (including dnSNVs and dnIndels) in 47 monozygotic twin pairs to validate the detection of EPZMs. Among 38 monozygotic pairs with available chorionicity and amnionicity data, 18.42% (7/38) were dichorionic diamniotic (DCDA), whereas 81.58% (31/38) were monochorionic diamniotic (MCDA). Across all 47 monozygotic pairs, we identified 65 shared EPZMs, along with 313 twin I-specific and 289 twin II-specific variants (Supplementary Table 18 ). This distribution indicates that 472 of the 602 EPZMs (78.41%) detected in individual twins likely arose after twinning. Comparative analysis revealed significantly higher EPZM concordance in MCDA twins than in DCDA twins ( P  = 0.019), supporting the hypothesis that later embryonic division (as in MCDA twins) leads to greater sharing of post-zygotic mutations. By contrast, analysis of 430 dizygotic pairs showed that only two shared EPZMs compared to 2,441 twin I-specific and 2,274 twin II-specific EPZMs (Supplementary Table 18 ). We employed a haplotype assembly approach to infer the parental origin of the DNM allele that occurred during gametogenesis 4 , 7 , 15 , 20 , 63 . In brief, considering that human autosomes are diploid and a DNM affects only one allele inherited from either the father or the mother, the other informative variants (uniquely attributable to a single parent) on the same haplotype could help distinguish the origin of the allele of the DNM. In addition, we applied a pedigree-based haplotype inference method 64 to increase the inference rate of informative variants. The inference rate is similar to those of phased mutations reported in recently published studies 4 , 7 , 15 , 20 , 63 . All dnSNVs were classified into three categories: (1) functional coding variants, including loss-of-function (LoF) mutations and missense variants predicted to be pathogenic by AlphaMissense; (2) gene-proximal variants, including non-pathogenic missense, synonymous, intronic and other variants of nearby genes; and (3) intergenic variants. To assess whether dnSNVs were enriched in regulatory elements, regulatory annotations were obtained from the Roadmap Epigenomics and ENCODE projects. Annotated regions included DNase I hypersensitive sites and histone modification marks associated with active regulatory elements. Analyses were performed across a diverse set of tissues and cell types, including ovary and testis tissues, six embryonic stem cell or trophoblast-derived cell lines (H1, H7, H9, HUES6, HUES48 and HUES64) and four differentiated embryonic germ-layer cell types. For baseline characteristics comparison across different groups, the Studentʼs t -test was used for continuous variables, and the χ 2 test was used for categorical variables. When accessing associations of DNMs with parental age, infant sex, ART, parental infertility diagnosis and ART-related procedures, DNMs were analyzed as standardized rates, obtained by normalizing raw DNM counts by callable genome size and scaling to the full genome length. These standardized DNM measures were modeled using generalized linear models (GLMs) with an identity link under a Gaussian assumption among singleton births. Given the non-independence of observations from the twin pairs, linear mixed-effects models (LMMs) with a family-level random intercept were additionally fitted to account for non-independence within twin pairs. Restricted cubic spline models were performed to assess the nonlinear associations of maternal DNMs with maternal age at conception. For analyses of offspring outcomes, GLMs were fitted for continuous outcomes (that is, gestational duration and birth weight), and logistic regression models were used for binary outcomes (that is, preterm birth, low birth weight and birth defects). All models adjusted for parental age at conception (continuous variable), maternal pre-pregnancy body mass index (BMI) (continuous variable), maternal education status, household income, area of residence, infant sex and mode of conception (ART or NC), with additional adjustment for feeding pattern and gestational duration (continuous variable) in analyses of neurodevelopmental outcomes. Model assumptions were evaluated using residual diagnostics and inspection of fitted versus observed values. Because only circulatory system birth defects (that is, CHD) had a prevalence greater than 1%, analyses of birth defects were restricted to this subtype (Extended Data Table 3 ). Mediation effects were estimated using a counterfactual-based framework implemented in the ‘mediation’ R package (version 4.5.0), with nonparametric bootstrap resampling with 1,000 iterations to derive 95% CIs for the mediated proportion. Continuous outcomes (that is, gestational duration and birth weight) were modeled using linear regression models, and binary outcomes (that is, preterm birth and low birth weight) were analyzed using logistic regression models. All mediation models were adjusted for the same set of covariates as the corresponding total-effect models. In all regression analyses, the outcome variable was the standardized absolute number of DNMs per infant, and, thus, all reported β coefficients represent differences in the absolute number of DNMs rather than standardized effect sizes. For continuous covariates (for example, parental age at conception), β corresponds to the expected change in the number of DNMs per 1-year increase in age. For binary covariates, β represents the mean difference in the number of DNMs between groups. No covariates were standardized to z -scores unless explicitly stated. Analyses were performed using R software (version 3.6.1), and a two-sided P  < 0.05 was considered to be statistically significant. Multiple testing correction was performed using Bonferroni correction where necessary. To minimize potential confounding arising from differences in sequencing platform, we repeated the primary ART−DNM association analyses restricted to samples sequenced exclusively on the NovaSeq platform (Supplementary Table 19 ). To further account for potential confounding by sequencing depth, we performed 1:2 matching between ART-conceived and NC infants based on parental and offspring sequencing depth using the R package MatchIt, and primary association models were then conducted in the depth-matched subset (Supplementary Table 20 ). The consistency of results across these platform-restricted and depth-matched sensitivity analyses indicates that the observed association between ART and increased DNM burden is unlikely to be driven by sequencing platform differences or coverage-related artifacts. Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

We performed WGS of 8,328 live-born infants from 7,851 completed parent−offspring families after quality control (Fig. 2 ), including 7,374 singleton and 477 twin pregnancies, yielding 24,030 analyzed samples. Among these families, 5,619 conceived naturally (NC) and 2,232 conceived through ART (Extended Data Table 1 ). Fig. 2 Flow diagram of our family-based WGS study. Flow diagram illustrating participant recruitment, quality-control procedures, exclusions and the final cohort included in the DNM analysis. QC, quality control. Source data Flow diagram illustrating participant recruitment, quality-control procedures, exclusions and the final cohort included in the DNM analysis. QC, quality control. Source data We identified 390,924 high-quality autosomal dnSNVs and 28,032 de novo indels (dnIndels), corresponding to an average mutation rate of 1.14 × 10 −8 per nucleotide per generation, within the estimated ranges between 1.00 × 10 −8 and 1.29 × 10 −8 from previous WGS studies 4 , 7 , 15 , 18 . Of these, 42,811 (10.22%; 39,330 dnSNVs and 3,481 dnIndels) with variant allele frequencies (VAFs) below 0.35 were classified as EPZMs. For the remaining 376,145 DNMs, parent of origin was determined for 39.18% DNMs (147,372), yielding 112,037 (29.79%; 104,954 dnSNVs and 7,083 dnIndels) paternal and 35,335 (9.39%; 33,040 dnSNVs and 2,295 dnIndels) maternal DNMs. In subsequent analyses, DNM counts were standardized by callable region size, resulting in an average standardized DNM count of 61.61 per offspring. Consistent with previous studies 4 , 19 , 20 , maternal dnSNVs showed enrichment of C > T transitions at non-CpG sites, whereas paternal dnSNVs were enriched for T > G, C > A, C > G and C > T substitutions at CpG sites ( P  < 3.57 × 10 −3 (0.05 / (7 × 2))) (Fig. 1b ). Trinucleotide context analysis further revealed diverse mutational spectra between paternal and maternal dnSNVs: maternal dnSNVs were significantly enriched for GTC > GAC and GTT > GAT, whereas paternal dnSNVs displayed preferential acquisition of ATG > AAG ( P   A and C > A substitutions and depletion of T > C substitutions (Fig. 1b ). Context-specific analysis further revealed that C > G transversions, especially GCG > GGG, CCG > CGG, TCG > TGG and ACG > AGG, were enriched in EPZMs ( P  < 2.60 × 10 −4 (0.05 / (96 × 2))) (Fig. 1c ). We then categorized dnSNVs into three classes. The vast majority of dnSNVs were noncoding, with 63.97% (250,075 / 390,924) as gene-proximal and 35.69% (139,532 / 390,924) as intergenic, whereas only 0.34% (1,317 / 390,924) were functional coding variants. These proportions were highly similar between paternal and maternal dnSNVs and between EPZMs and germline dnSNVs (Fisherʼs exact test, P  > 8.33 × 10 −3 (0.05 / (3 × 2)); Fig. 1d ), indicating no enrichment of functional variants by mutational origin. Similarly, no enrichment was observed across regulatory annotations, including histone-marked regions and DNase I hypersensitive sites (Supplementary Fig. 1 ). Finally, we assessed the clinical relevance and pathogenicity of dnSNVs according to the American College of Medical Genetics and Genomics (ACMG) V3 (ref. 21 ) draft V4 and current Association for Clinical Genome Sequencing (ACGS) guidelines 22 . Across the entire dataset (Supplementary Table 1 ), only three dnSNVs were classified as pathogenic or likely pathogenic (P/LP) by ClinGen, including confirmed RASopathy-associated mutations in SOS1 (p.Arg552Ser) 23 and PTPN11 (p.Met504Val) 24 and a PAH (p.Pro407Leu) variant 25 . An additional three dnSNVs were classified as high-confidence P/LP (at least a two-star rating) within ClinVar. Although rare, their identification highlights the pathogenic potential of de novo mutational process. Consistent with the ‘male bias’ observed in previous studies 4 – 8 , we observed a stronger paternal than maternal age effect after accounting for both parental ages at conception in 7,374 singletons (Fig. 3a ). Each additional year of paternal age was associated with an increase of 1.06 paternal DNMs (95% confidence interval (CI): 0.99−1.14, P  = 2.15 × 10 −164 ), including 1.00 dnSNVs (95% CI: 0.93−1.07, P  = 1.68 × 10 −159 ) and 0.058 dnIndels (95% CI: 0.038−0.078, P  = 7.80 × 10 −9 ). By contrast, each additional year of maternal age was associated with 0.38 maternal DNMs (95% CI: 0.32−0.44, P  = 2.11 × 10 −33 ), including 0.36 dnSNVs (95% CI: 0.30−0.42, P  = 1.27 × 10 −32 ) and 0.017 dnIndels (95% CI: 0.003−0.031, P  = 0.021). Fig. 3 The effect of parental age at conception on DNMs. a , Associations between parental age at conception and the standardized number of parental DNMs among 7,374 singleton births. Lines represent fitted estimates from GLMs including parental age at conception, shaded areas indicate the 95% CIs and regression coefficients ( β ) represent the estimated change in the number of DNMs per 1-year increase in parental age. b , Absolute parental age effects for each mutational class among 7,374 singleton births. Points indicate regression coefficients ( β ), and error bars indicate the 95% CIs. The line represents a linear regression through the origin using the numbers of paternal DNMs for the mutational classes as weights. c , Restricted cubic spline plots of the association between maternal age at conception and the standardized number of maternal dnSNVs among 7,374 singleton births, adjusted for paternal age at conception. The solid line indicates the fitted estimate; the shaded area indicates the 95% CI; and the histogram shows the distribution of maternal age at conception. d , Relative frequency of C > G substitutions at the APOBEC-related TCA or TCT trinucleotide contexts by maternal age at conception. Proportions were compared using two-sided Fisherʼs exact tests. ** P  < 0.01; *** P  < 0.001. Exact P values are provided in the Source Data file. Source data a , Associations between parental age at conception and the standardized number of parental DNMs among 7,374 singleton births. Lines represent fitted estimates from GLMs including parental age at conception, shaded areas indicate the 95% CIs and regression coefficients ( β ) represent the estimated change in the number of DNMs per 1-year increase in parental age. b , Absolute parental age effects for each mutational class among 7,374 singleton births. Points indicate regression coefficients ( β ), and error bars indicate the 95% CIs. The line represents a linear regression through the origin using the numbers of paternal DNMs for the mutational classes as weights. c , Restricted cubic spline plots of the association between maternal age at conception and the standardized number of maternal dnSNVs among 7,374 singleton births, adjusted for paternal age at conception. The solid line indicates the fitted estimate; the shaded area indicates the 95% CI; and the histogram shows the distribution of maternal age at conception. d , Relative frequency of C > G substitutions at the APOBEC-related TCA or TCT trinucleotide contexts by maternal age at conception. Proportions were compared using two-sided Fisherʼs exact tests. ** P  < 0.01; *** P  < 0.001. Exact P values are provided in the Source Data file. Source data Analysis by mutation classes revealed substantial heterogeneity in the maternal age effect of dnSNVs (Fig. 3b and Supplementary Table 2 ). The paternal age effect could be observed for all mutational classes ( P   T transitions at CpG sites ( P  = 0.30), consistent with previous evidence that deamination-derived CpG C > T mutations accumulate less extensively in the female germline 26 . Despite differences in mutational spectra, the two most common types affected by parental age were T > C and C > T substitutions at non-CpG sites. In addition, although the ‘male bias’ was stable with parental age for most mutation types, maternal age effects on C > G and T > A were similar to those of paternal dnSNVs. Results were consistent when twin pregnancies were included (Supplementary Table 3 ). Then, we considered polynomial age effects and observed that a quadratic model provided a better fit than the linear model for maternal, but not paternal, dnSNVs (maternal age: difference in the Akaike information criterion (ΔAIC) = 8.44, P  = 1.23 × 10 −3 ; paternal age: ΔAIC = 1.83, P  = 0.68) (Supplementary Table 4 ). Then, we performed restricted cubic spline models with three knots at the 50th (approximately 28.92 years), 60th (approximately 29.94 years) and 75th (approximately 31.81 years) percentiles of maternal age, and the estimated slope in mothers over 31.81 years of age (the old group) was approximately 1.52-fold of that of mothers under 28.92 years of age (the young group) after accounting for paternal age (0.42 versus 0.28, P for nonlinear test = 1.05 × 10 −3 ) (Fig. 3c ). This pattern remained after excluding offspring with zero maternally assignable dnSNVs and after rank-based inverse normal transformation of maternal dnSNVs (Supplementary Fig. 2 ). In addition, no evidence of nonlinearity was observed in the young ( P for nonlinear test = 0.61) and old ( P for nonlinear test = 0.70) groups, suggesting an accelerated maternal age effect around the age between 28.92 and 31.81. To explore the mechanisms underlying the quadratic maternal age effect, we compared linear and quadratic models across seven mutational classes. A quadratic maternal age effect provided a significantly better fit only for C > G transversions ( P for nonlinear test = 7.45 × 10 −3 ), whereas linear models adequately described the remaining six mutation types ( P for nonlinear test > 0.05) (Extended Data Fig. 1 ). Further analyses of 16 trinucleotide contexts flanking C > G substitutions identified that the two most significantly enriched mutation contexts in the old maternal age group (>31.81 years) were T(C > G)T (Fisherʼs exact test P  = 2.50×10 −4 ) and T(C > G)A (Fisherʼs exact test P  = 6.30×10 −3 ) (Fig. 3d and Supplementary Table 5 ). The enrichment of TCA/TCT trinucleotide contexts remained significant, even after adjustment for mutational class composition 4 ( P  = 0.040; Supplementary Table 6 ). Although maternal C > G mutations have previously been attributed to non-crossover gene conversions in aging oocytes 4 , 27 , the enrichment of canonical APOBEC-associated TCA/TCT signatures suggests alternative mechanisms in the germline of older mothers, such as REV1-dependent translesion synthesis (TLS) 28 , 29 . In addition to parental age at conception, we investigated associations between ART and the standardized number of DNMs. After adjusting for parental age at conception and infant sex, ART-conceived singleton infants carried 2.48 more DNMs than NC infants (95% CI: 1.94−3.02, P  = 2.53 × 10 −19 ), and this association was primarily driven by dnSNVs rather than dnIndels (Extended Data Table 2 ). Specifically, ART-conceived singleton infants carried 2.38 more dnSNVs than NC infants (95% CI: 1.87−2.89, P  = 8.73 × 10 −20 ), including 0.74 paternal dnSNVs (95% CI: 0.25−1.22, P  = 2.96 × 10 −3 ), 0.62 maternal dnSNVs (95% CI: 0.28−0.96, P  = 3.30 × 10 −4 ) and 1.02 EPZMs (95% CI: 0.86−1.18, P  = 7.03 × 10 −34 ), indicating that nearly half (1.02/2.38 = 42.86%) of ART-associated dnSNVs occur during early embryonic development (Fig. 4a ). Fig. 4 Association between detailed ART procedures and germline dnSNVs. a , ART-conceived infants ( n  = 1,782) showed significantly more dnSNVs than NC infants ( n  = 5,592), with adjustment for parental age at conception and infant sex. Box plots display the median, IQR and whiskers within the range of 1.5× IQR. Regression coefficients ( β ) were estimated using a GLM including parental age at conception and infant sex. b , Associations between ART procedural parameters and dnSNV burden estimated using a GLM including parental age at conception, infant sex, fertilization method, embryo transfer stage, embryo transfer cycle, ovulation induction regimen and trigger strategy. c , Association between fertilization method and paternal dnSNV burden estimated using a GLM adjusted for parental age at conception and infant sex. d , Association between the dosage of the regimen of GnRH-ant and maternal dnSNV burden estimated using a GLM adjusted for parental age at conception and infant sex. For b − d , squares indicate regression coefficients ( β ), and error bars indicate the 95% CIs. Exact sample sizes for each ART procedure category are provided in the Source Data file. GnRHa, gonadotropin-releasing hormone agonist. Source data a , ART-conceived infants ( n  = 1,782) showed significantly more dnSNVs than NC infants ( n  = 5,592), with adjustment for parental age at conception and infant sex. Box plots display the median, IQR and whiskers within the range of 1.5× IQR. Regression coefficients ( β ) were estimated using a GLM including parental age at conception and infant sex. b , Associations between ART procedural parameters and dnSNV burden estimated using a GLM including parental age at conception, infant sex, fertilization method, embryo transfer stage, embryo transfer cycle, ovulation induction regimen and trigger strategy. c , Association between fertilization method and paternal dnSNV burden estimated using a GLM adjusted for parental age at conception and infant sex. d , Association between the dosage of the regimen of GnRH-ant and maternal dnSNV burden estimated using a GLM adjusted for parental age at conception and infant sex. For b − d , squares indicate regression coefficients ( β ), and error bars indicate the 95% CIs. Exact sample sizes for each ART procedure category are provided in the Source Data file. GnRHa, gonadotropin-releasing hormone agonist. Source data Mutation-class analyses showed ART-associated effects across all EPZM mutation classes, whereas parent-of-origin dnSNVs exhibited more specific patterns (Supplementary Figs. 3 and 4 ). Additional analyses of pathogenic variants in core DNA damage repair genes showed no enrichment among ART parents, and adjustment for parental carrier status did not materially change the observed ART-associated effects (Supplementary Tables 7 – 9 ). These results do not support the hypothesis that elevated baseline germline mutational susceptibility due to impaired DNA repair capacity explains the observed differences in offspring dnSNV burden. We then explored whether specific parental infertility diagnoses or ART procedures were associated with dnSNVs among ART-conceived births after adjustment for parental age at conception and infant sex (Supplementary Tables 10 and 11 ). ART procedures showing evidence of association with dnSNV burden, including gonadotropin-releasing hormone antagonist (GnRH-ant) use, trigger strategy, frozen embryo transfer, blastocyst-stage transfer and ICSI, were carried forward into a prespecified multivariable model. By contrast, no consistent evidence of association was observed for infertility diagnoses, including oligoasthenoteratozoospermia (OAT), endometriosis, ovarian dysfunction, polycystic ovary syndrome or infertility due to pelvic or tubal factor ( P  > 0.05). Similar findings were observed when infertility diagnoses were aggregated into broader composite categories (Supplementary Table 12 ). We next fitted a multivariable regression model in which the ART procedural parameters identified in the preliminary screening analyses were further included together with parental age at conception and infant sex. In this model, ART procedures at the gamete stage were specifically associated with increased parent-of-origin dnSNVs. Specifically, ICSI was associated with elevated paternal dnSNVs ( β  = 1.01, 95% CI: 0.02−2.00, P  = 0.046), whereas GnRH-ant use was associated with increased maternal dnSNVs ( β  = 1.69, 95% CI: 0.49−2.89, P  = 5.95 × 10 −3 ) (Fig. 4b and Supplementary Table 13 ). These associations remain consistent when twin pregnancies were included (Supplementary Fig. 5a and Supplementary Table 14 ). As ICSI is routinely employed for couples with paternal OAT, OAT was substantially more common among ICSI than in vitro fertilization (IVF) conceptions (70.2% (318/453) versus 39.1% (479/1,224); Fisherʼs exact test P  = 5.05 × 10 −30 ), and all 44 men diagnosed with azoospermia underwent ICSI. We conducted additional stratified analyses based on paternal OAT status in ART-conceived infants. ICSI was consistently associated with increased paternal dnSNVs than IVF, regardless of OAT diagnosis ( P for heterogeneity = 0.83) (Supplementary Fig. 6 ). When the NC infants were further included (Fig. 4c ), paternal dnSNVs in IVF conceptions were similar to the NC group ( β  = 0.41, 95% CI: −0.15 to 0.97, P  = 0.15), whereas ICSI conceptions exhibited significantly higher dnSNVs than the NC group ( β  = 1.58, 95% CI: 0.74–2.43, P  = 2.54 × 10 −4 ), suggesting that ICSI largely accounted for the excess paternal DNMs observed among ART-conceived infants. We then evaluated the impact of GnRH-ant dosage on maternal dnSNVs and stratified families with the usage of GnRH-ant into three dosing groups based on the current standard regimen of 0.25 mg per day: the half-dose group (≤0.125 mg per day), the full-dose group (0.125−0.25 mg per day) and the high-dose group (>0.25 mg per day). A significant dose−response correlation was observed ( P for trend = 9.41 × 10 −3 ): for women in the high-dose group, their children carried more dnSNVs than those who did not use the GnRH-ant protocol (Supplementary Fig. 7 ). When the NC group was further included (Fig. 4d ), maternal dnSNVs among ART conceptions without GnRH-ant treatment were similar to the NC group ( β  = 0.21, 95% CI: −0.25 to 0.66, P  = 0.37), whereas the maternal dnSNVs among the GnRH-ant-treated group were significantly higher, particularly in the high-dose group ( β  = 1.65, 95% CI: 0.08−3.23, P  = 0.040), indicating that GnRH-ant treatment may contribute to the increase of maternal dnSNVs in ART-conceived infants. Results were consistent when twin pregnancies were included (Supplementary Fig. 5b,c ). We next evaluated the contribution of parental dnSNVs to parental age-related and ART-related offspring health outcomes among 7,374 singleton births (Fig. 5a and Extended Data Table 3 ). After adjusting for parental ages and other potential confounders ( Methods ), each additional paternal dnSNV was associated with shorter gestational duration (week) ( β  = −0.006, 95% CI: −0.010 to −0.002, P  = 1.78 × 10 − 3 ) and lower birth weight (per 100 g; β  = −0.018, 95% CI: −0.029 to −0.006, P  = 2.68 × 10 − 3 ) as well as increased risks of preterm birth (odds ratio = 1.02, 95% CI: 1.01−1.04, P  = 2.39 × 10 − 4 ) and low birth weight (odds ratio = 1.02, 95% CI: 1.01−1.04, P  = 3.27 × 10 − 3 ) ( P   0.05). Notably, a significant association was detected only for functional dnSNVs in established congenital heart disease (CHD) genes (Supplementary Tables 15 and 16 and Supplementary Fig. 8 ), which were associated with a 39-fold increased risk of CHD (odds ratio = 39, 95% CI: 6−246, P  = 9.30 × 10 −5 ) (Fig. 5b ). Fig. 5 Association between parent-of-origin dnSNVs and offspring health outcomes. a , Association between paternal and maternal dnSNVs and offspring health outcomes among 7,374 singleton births. Regression coefficients ( β ) for continuous outcomes were estimated using GLMs, and odds ratios (ORs) for binary outcomes were estimated using logistic regression models, adjusted for parental ages and other potential confounders described in the Methods section. For continuous outcomes, summary values are shown as mean ± s.d.; for binary outcomes, prevalence is shown as n (%). Squares indicate β coefficients or ORs, and error bars indicate the 95% CIs. b , Association of the overall burden of parental dnSNVs (dnSNVs P ) and functional parental dnSNVs (fdnSNVs P ) with risk of CHD among 7,374 singleton births, estimated using logistic regression models adjusted for parental age at conception, maternal pre-pregnancy BMI, maternal education status, household income, area of residence, infant sex and mode of conception. Points indicate ORs, and error bars indicate the 95% CIs. fdnSNVs P were defined as LoF variants and P/LP variants defined in the ClinVar dataset. Exact P values are provided in the Source Data file. c , Mediation analysis assessing the mediating role of paternal dnSNVs in the association between parental age and gestational duration. d , Mediation analysis assessing the mediating role of paternal dnSNVs in the association between ICSI and gestational duration. For c and d , mediated proportions were estimated using the ‘mediation’ R package with 1,000 nonparametric bootstrap iterations to derive 95% CIs. Solid lines indicate statistically significant paths. * P  < 0.05; ** P  < 0.01; *** P  < 0.001. Exact P values are provided in the Source Data file. Prop., proportion. Source data a , Association between paternal and maternal dnSNVs and offspring health outcomes among 7,374 singleton births. Regression coefficients ( β ) for continuous outcomes were estimated using GLMs, and odds ratios (ORs) for binary outcomes were estimated using logistic regression models, adjusted for parental ages and other potential confounders described in the Methods section. For continuous outcomes, summary values are shown as mean ± s.d.; for binary outcomes, prevalence is shown as n (%). Squares indicate β coefficients or ORs, and error bars indicate the 95% CIs. b , Association of the overall burden of parental dnSNVs (dnSNVs P ) and functional parental dnSNVs (fdnSNVs P ) with risk of CHD among 7,374 singleton births, estimated using logistic regression models adjusted for parental age at conception, maternal pre-pregnancy BMI, maternal education status, household income, area of residence, infant sex and mode of conception. Points indicate ORs, and error bars indicate the 95% CIs. fdnSNVs P were defined as LoF variants and P/LP variants defined in the ClinVar dataset. Exact P values are provided in the Source Data file. c , Mediation analysis assessing the mediating role of paternal dnSNVs in the association between parental age and gestational duration. d , Mediation analysis assessing the mediating role of paternal dnSNVs in the association between ICSI and gestational duration. For c and d , mediated proportions were estimated using the ‘mediation’ R package with 1,000 nonparametric bootstrap iterations to derive 95% CIs. Solid lines indicate statistically significant paths. * P  < 0.05; ** P  < 0.01; *** P  < 0.001. Exact P values are provided in the Source Data file. Prop., proportion. Source data We then evaluated the contribution of parental dnSNVs to parental age effects specifically on birth outcomes. Because maternal and paternal ages were highly correlated ( r  = 0.78, P  < 0.001), we adopted two complementary modeling strategies. In one model, we used the mean parental age to represent the shared age-related component of both parents. In a separate model, we simultaneously included maternal and paternal ages. Average parental age was significantly associated with shorter gestational duration and increased risk of low birth weight (Extended Data Fig. 2 ). Mediation analyses, therefore, focused on these two outcomes and indicated that paternal dnSNVs statistically accounted for 13.57% (95% CI: 3.4−25%, P  = 0.01) of the parental age effect on gestational duration shortening (Fig. 5c ) and for 44.70% (95% CI: 14.6−147%, P  = 0.004) of the effect on low birth weight (Extended Data Fig. 3 ) ( P  < 0.025 (0.05 / 2)). No significant mediation effects were expected for the independent paternal age or maternal age effects (Fig. 5a and Extended Data Fig. 4 ). Although average parental age was also associated with several neurodevelopmental delays (Extended Data Fig. 2 ), none of these was associated with parent-of-origin dnSNVs (Fig. 5a ), precluding further mediation analyses. We then examined ART-related offspring outcomes (Extended Data Fig. 5 ). Compared to NC infants, ART-conceived infants had significantly shorter gestational duration ( β  = −0.33, 95% CI: −0.41 to −0.25, P  = 8.42 × 10 − 16 ) and lower birth weight ( β  = −0.44, 95% CI: −0.69 to −0.19, P  = 6.49 × 10 − 4 ) as well as higher risks of preterm birth (odds ratio = 1.81, 95% CI: 1.41−2.33, P  = 3.51 × 10 − 6 ) and low birth weight (odds ratio = 1.62, 95% CI: 1.18−2.23, P  = 2.74 × 10 − 3 ) after accounting for potential confounders ( Methods ). Because ART was also associated with increased paternal dnSNVs, mediation analyses were performed to quantify the contribution of ART-associated dnSNVs to birth outcomes. Paternal dnSNVs explained 1.29% (95% CI: 0.23−3%, P  = 0.008) of the ART-associated reduction in gestational duration and 2.95% (95% CI: 0.59−7%, P  = 0.008) of the excess risk of preterm birth ( P  < 0.0125 (0.05 / 4)) (Extended Data Fig. 6a−c ). Because the observed increase in paternal dnSNVs in ART conceptions was primarily attributable to ICSI, we further compared NC and ICSI births, in which paternal dnSNVs mediated 3.99% (95% CI: 0.77−14%, P  = 0.014) of the difference in gestational duration (Fig. 5d ) and 8.32% (95% CI: 1.44−42%, P  = 0.034) of the excess preterm birth risk (Extended Data Fig. 6d ). Although ART-conceived infants also exhibited higher risks of CHD and delays in gross motor and cognitive development (Extended Data Fig. 5 ), these outcomes were not associated with parent-of-origin dnSNVs, precluding further mediation analyses. Unlike parental dnSNVs, EPZMs were specifically associated with embryo transfer cycle and transfer stage (Supplementary Tables 10 and 11 ). In the multivariate model including other ART procedures, parental age at conception and infant sex, frozen embryo transfer was associated with an increased number of EPZMs in ART-conceived singletons ( β  = 0.61, 95% CI: 0.13−1.09, P  = 0.012) (Supplementary Table 13 ). Analyses including twin pregnancies revealed that both embryo transfer cycle and transfer stage had associations with EPZMs (blastocyst versus cleavage stage: β  = 0.31, 95% CI: 0.01−0.62, P  = 0.045; frozen versus fresh embryo: β  = 0.37, 95% CI: −0.01 to 0.76, P  = 0.054) (Supplementary Table 14 ). We subsequently categorized the study population into four groups based on transfer cycle type and transfer stage. Compared to fresh cleavage-stage embryo transfer, frozen blastocyst-stage embryo transfer was associated with 0.81 more EPZMs (95% CI: 0.31−1.30, P  = 1.37 × 10 −3 ), where fresh blastocyst-stage embryo transfer and frozen cleavage-stage embryo transfer showed intermediate increases (Supplementary Fig. 9a,b ). When the NC group was included (Fig. 6a ), both transfer cycle and transfer stage were associated with higher EPZM burdens than the NC group. In particular, frozen blastocyst-stage embryo transfer was associated with 1.23 more EPZMs compared to NC infants (95% CI: 1.02−1.44, P  = 1.38 × 10 −30 ). Results were consistent when twin pregnancies were included (Supplementary Fig. 9c ). These findings suggest that early embryonic manipulations during ART may influence the generation of EPZMs. Fig. 6 Association between ART and EPZMs. a , Association of embryo transfer stage and transfer cycle with EPZM burden estimated using GLMs, adjusted for parental age at conception and infant sex. Squares indicate regression coefficients ( β ), and error bars indicate the 95% CIs. The right panel shows the proportion of C > A substitutions among EPZMs for each group, with percentages labeled beside each bar. Exact P values are provided in the Source Data file. b , Association of specific mutation classes with embryo transfer stage and embryo transfer cycle. Estimates were derived from GLMs with adjustment for parental age at conception and infant sex. Squares indicate regression coefficients ( β ), and error bars indicate the 95% CIs. Exact sample sizes for each ART procedure category are provided in the Source Data file. c , Comparative analyses of proportions of the mutational classes of EPZMs among the NC group and the ART group stratified by embryo transfer stage or embryo transfer cycle. Proportions were compared using two-sided Fisherʼs exact tests. Exact P values are provided in the Source Data file. d , Correlation of the differential mutational profile between blastocyst-transfer and NC offspring with COSMIC SBS signatures. SBS18, a signature associated with ROS-induced oxidative DNA damage, showed the strongest similarity to the observed mutational profile. Exact correlation coefficients and P values are provided in the Source Data file. e , Association between C > A EPZMs and the risk of delays in neurodevelopment at 1 year of age among 7,374 singleton births. Estimates were derived from logistic regression models with adjustment for parental age at conception, maternal pre-pregnancy BMI, maternal education status, household income, area of residence, infant sex, mode of conception, feeding pattern and gestational duration. Squares indicate odds ratios (ORs), and error bars indicate the 95% CIs. For a , c and d , * P  < 0.05; ** P  < 0.01; *** P  < 0.001. n , number of infants. Source data a , Association of embryo transfer stage and transfer cycle with EPZM burden estimated using GLMs, adjusted for parental age at conception and infant sex. Squares indicate regression coefficients ( β ), and error bars indicate the 95% CIs. The right panel shows the proportion of C > A substitutions among EPZMs for each group, with percentages labeled beside each bar. Exact P values are provided in the Source Data file. b , Association of specific mutation classes with embryo transfer stage and embryo transfer cycle. Estimates were derived from GLMs with adjustment for parental age at conception and infant sex. Squares indicate regression coefficients ( β ), and error bars indicate the 95% CIs. Exact sample sizes for each ART procedure category are provided in the Source Data file. c , Comparative analyses of proportions of the mutational classes of EPZMs among the NC group and the ART group stratified by embryo transfer stage or embryo transfer cycle. Proportions were compared using two-sided Fisherʼs exact tests. Exact P values are provided in the Source Data file. d , Correlation of the differential mutational profile between blastocyst-transfer and NC offspring with COSMIC SBS signatures. SBS18, a signature associated with ROS-induced oxidative DNA damage, showed the strongest similarity to the observed mutational profile. Exact correlation coefficients and P values are provided in the Source Data file. e , Association between C > A EPZMs and the risk of delays in neurodevelopment at 1 year of age among 7,374 singleton births. Estimates were derived from logistic regression models with adjustment for parental age at conception, maternal pre-pregnancy BMI, maternal education status, household income, area of residence, infant sex, mode of conception, feeding pattern and gestational duration. Squares indicate odds ratios (ORs), and error bars indicate the 95% CIs. For a , c and d , * P  < 0.05; ** P  < 0.01; *** P   A EPZMs exhibited a significant association with blastocyst-stage embryo transfer ( β  = 0.13, 95% CI: 0.04−0.22, P  = 3.86 × 10 −3 ) after adjusting for parental age and infant sex ( P   A EPZMs was significantly higher in blastocyst-stage compared to cleavage-stage transfers (13.45% versus 12.28%, P  = 0.047) and natural conceptions (13.45% versus 11.95%, P  = 5.24 × 10 −4 ), whereas cleavage-stage transfers showed similar rates to natural conceptions ( P  = 0.48) (Fig. 6c ). Similarly, frozen embryo transfers exhibited a higher proportion of C > A EPZMs than natural conceptions (13.20% versus 11.95%, P  = 8.98 × 10 −4 ), whereas fresh embryo transfers did not (12.04% versus 11.95%, P  = 0.90) (Fig. 6c ). Because C > A transversions are mechanistically associated with 8-oxoguanine (8-oxoG) and 8-oxodeoxyguanosine 30 , common DNA lesions caused by reactive oxygen species (ROS), we propose that the increased C > A transversions in blastocysts may stem from oxidative DNA damage during prolonged culture. Supporting this hypothesis, mutational signatures distinguishing blastocyst transfers from natural conceptions showed the strongest similarity to COSMIC single base substitution (SBS) signature 18 (SBS18) ( P  = 8.75 × 10 − 3 ) (Fig. 6d ), a signature linked to ROS-induced DNA damage 31 . Given the well-known association between post-zygotic mutations and neurodevelopment, we evaluated the association between C > A EPZMs and delayed neurodevelopment. Each additional C > A EPZM was associated with a 14% higher risk of delayed neurocognitive development (95% CI: 1.01−1.28, P  = 0.035) (Fig. 6e ), whereas no significant association was observed for the overall burden of EPZM (Supplementary Fig. 10 ).

Through WGS of this large prospective birth cohort, we systematically evaluated whether parental age and ART shape mutational burden and whether DNM burden relates to offspring outcomes. We found that ART contributes to DNM accumulation independent of parental age at conception and identified stage-specific associations, where ICSI, ovarian stimulation and embryo culture were associated with paternal, maternal and early post-zygotic mosaic DNMs, respectively. Furthermore, by integrating longitudinal follow-up data, we demonstrated that paternal DNMs statistically accounted for the effects of both ICSI and advanced parental age on shortened gestational duration, whereas embryo culture-associated EPZMs were associated with delayed neurocognitive development. These findings provide additional insights into the origins and consequences of DNMs. Consistent with previous studies 7 , 20 , paternal age was the dominant determinant of DNM burden. Although DNMs have long been proposed as a molecular basis for the paternal age effect on offspring health 32 , genetic modeling studies suggest that age-related DNMs alone are unlikely to fully explain this effect 33 . Our findings provide additional population-based epidemiological evidence that paternal age does not influence offspring health directly via DNMs but, rather, through DNM-associated early-life traits, such as gestational duration and birth weight 34 , thereby offering an alternative explanation for the broader paternal age-related effects observed across diverse phenotypes. However, we did not observe mediation for congenital defects or neurodevelopmental traits, likely reflecting the general population nature of our cohort; this merits further evaluation in disease-specific populations 35 . Beyond point mutations, other genomic mechanisms 36 , including de novo structural variants and copy number alterations, also warrant consideration, yet their rarity and severe effects likely favor contributions to spontaneous disorders over complex diseases. Maternal DNMs also exhibit an age-dependent increase but display distinct mutational patterns. The relative depletion of CpG > TpG substitutions, together with their weak age dependence observed in our study, suggests that time-dependent deamination of methylated cytosines 37 contributes less to maternal mutagenesis than to paternal mutagenesis. By contrast, the enrichment and age-accelerated increase of T C A / T C T substitutions among older women, which are known APOBEC-associated signatures linked to REV1-dependent TLS 29 , point to an age-related remodeling of DNA damage tolerance or repair fidelity during oogenesis. Together, these findings indicate that maternal mutagenesis is governed by mechanisms distinct from those operating in paternal gametogenesis and cannot be explained solely by cell division-related processes. ART now accounts for more than 5% of births worldwide 38 , yet its mutagenic effects and health consequences in humans remain largely uncertain. Leveraging a prospective cohort with detailed ART procedural information, our study indicates that ART should not be treated as a single exposure but, rather, as a collection of interventions with potentially different mutational consequences. ICSI is now the most commonly used fertilization method in ART 39 . The elevated paternal dnSNVs observed among ICSI-treated conceptions may result from bypassing natural selection at fertilization 32 , 40 , 41 , thereby permitting the transmission of sperm carrying ‘selfish mutations’, although mechanical or biochemical stress during the ICSI procedure may also contribute. Because paternal dnSNVs did not differ by male-factor infertility, infertility alone is unlikely to account for the observed association. Interestingly, similar to paternal age, ICSI-associated paternal dnSNVs partially mediated the reduction in gestational duration, implying broader consequences for offspring health. For maternal dnSNVs, we observed a specific dose-dependent association with GnRH-ant use. Although GnRH-ant regimens offer several clinical advantages 42 – 44 , some studies have reported lower pregnancy and live birth rates than GnRH-agonist protocols 43 , 45 , 46 . Our findings provide a potential genomic perspective on such adverse effects. GnRH-induced meiotic maturation or related cellular stress during oogenesis 47 , 48 may affect genomic instability and increase maternal dnSNVs, and the dose-dependent effect further supports a biological rather than a stochastic effect. These findings suggest that individualized low-dose GnRH-ant strategies warrant further evaluation for both reproductive and potential mutational consequences. Benchmarking ART-related germline DNM shifts against natural parental aging provides a metric to gauge the biological cost of reproductive interventions. Using the estimated age-related increase in dnSNVs as a reference, the mutational effects of ICSI and GnRH-ant use were equivalent to 1.01 and 4.69 additional years of paternal and maternal aging, respectively. Because these germline DNMs are constitutional and heritable, their potential to alter the genetic landscape over generations highlights the importance of considering the cumulative and incalculable long-term impact on public health. Beyond parent-of-origin DNMs, we found that embryo freezing and preimplantation culture may influence mutational processes during early embryogenesis. As extended embryo culture becomes increasingly common in ART 49 , concerns have emerged because the embryo culture window is highly susceptible to mutagenesis, especially for post-zygotic mutations 17 . Our analysis suggests that oxidative stress induced by prolonged in vitro embryo culture 50 may explain the elevated C > A transversions observed in blastocyst transfers. Although the precise biological mechanisms require further investigation, our findings highlight opportunities to optimize culture conditions and refine in vitro culture strategies. Although our prospective cohort design strengthens the investigation into the generation and health effects of DNMs, several limitations deserve consideration. First, our focus on EPZMs with VAF less than 0.35 might underestimate the full burden of post-zygotic mutations and cannot capture tissue-specific mosaicism. Second, the restriction to live-born infants precludes a comprehensive assessment of fetal-lethal DNMs in cases of stillbirth or pregnancy termination, potentially leading to underestimation of DNM-associated adverse effects. Third, as with all observational studies, residual confounding cannot be fully excluded despite extensive adjustment for parental, clinical and procedural factors. Therefore, our results should be interpreted as population-level associations that generate biological insight rather than causal estimates or guidance for individual clinical decision-making. In summary, our study delineates how biological aging and ART shape de novo mutagenesis and influence offspring health in a large prospective birth cohort. By resolving parent-of-origin DNMs and EPZMs, we demonstrate that distinct stages of ART are associated with different mutational processes and that these mutations contribute to variation in offspring outcomes. Collectively, these findings elucidate the contribution of parent-of-origin DNMs and EPZMs to offspring health and highlight the importance of optimizing ART procedures to mitigate their mutational and developmental consequences.

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41591-026-04585-2 .

Source Data Figs. 1−6, Source Data Extended Data Figs. 2−6 and Source Data Extended Data Table 2. Statistical source data. Statistical source data.