cfGWAS reveal genetic basis of cell-free DNA end motifs

Domain: neuroscience
Published: 2026-02-14
Editorial outcome: major_revision

Abstract
Cell-free DNA (cfDNA) consists of degraded DNA fragments released into body fluids. Its genetic and pathological information makes it useful for prenatal testing and early tumor detection. However, the mechanisms behind cfDNA biology are largely unknown. In this study, for the first time, we conduct a genome-wide association study (GWAS) to explore the genetic basis of cfDNA end motif frequencies, termed cfGWAS, in 28,016 pregnant women. We identify 15 study-wide significant loci, including the well-known cfDNA-related genes DFFB and DNASE1L3, as well as novel genes potentially involved in cfDNA biology, such as PANX1 and DNASE1L1. The findings are further verified through three independent GWAS studies and experimental validation in knockout mice and cell lines. Subsequent analyses reveal strong causal relationships of leukocytes, especially neutrophils, with cfDNA features. In summary, we introduce the cfGWAS, revealing the genetic basis of cfDNA biology on a genome-wide scale. Novel knowledge uncovered by this study promises to revolutionize liquid biopsy technology and lead to potential new drugs targeting certain diseases. Given that millions of cfDNA whole genome sequencing data have been generated from clinical testing, the potential of this paradigm is enormous.

Introduction
Cell-free DNA (cfDNA) is DNA fragments released from cells, often originating from processes such as apoptosis, necrosis, active secretion, and the specialized cell death process known as NETosis1,2. These DNA fragments circulate freely in bodily fluids such as blood, saliva, and urine, reflecting the biological status of cells or tissues. CfDNA analysis has extensive applications in non-invasive prenatal testing (NIPT) for detecting fetal abnormalities3. It has also emerged as a valuable biomarker for early cancer screening and monitoring4,5, as well as for detecting organ rejection in transplantation6. Abnormalities in cfDNA have been observed in autoimmune diseases such as systemic lupus erythematosus (SLE)7. Additionally, cfDNA shows promise as a tool for diagnosing and managing infectious diseases8,9.

Commonly used characteristics of cfDNA include cfDNA concentration, fragment length, end motifs, jagged ends, as well as nucleosome footprints10. In healthy individuals, cfDNA concentration is low, but it significantly increases in diseases like cancer, inflammation, and tissue damage11–14. Therefore, measuring cfDNA concentration can help with disease prediction and monitoring. The fragment length of plasma cfDNA fragments typically ranges from 50 to 600 base pairs, though both shorter and longer fragments have also been observed15,16, but this may vary due to individual differences, physiological states, and disease conditions10. Researchers analyzed the fragment lengths to capture the disease signals, such as aberrant shortening in cfDNA from cancer patients17,18. End motifs refer to the nucleotide composition at the 5’ end of cfDNA, generated during cfDNA digestion. Different cfDNA end motifs reflect nuclease activity, pathophysiological processes, and tissue-specific DNA fragmentation patterns. Specifically, DFFB, DNASE1L3, and DNASE1 preferentially generate A-end, C-end, and T-end motifs, respectively11. Recent studies have highlighted cfDNA end motifs as emerging biomarkers for oncology and transplantation monitoring due to their frequency shifts under pathological conditions and their tissue-specific profiles19.

These cell-free DNA features arise from the degradation of DNA both inside and outside of cells, and some may correlate due to shared biological processes10,20,21. Understanding the genes responsible for cell-free DNA fragmentation and their functions is essential for elucidating the biological mechanisms that influence cfDNA fragmentomic profiles in both healthy and disease states. Currently, there is a lack of effective methods for identifying these genes, which has restricted previous studies to examining those associated with cell death and nuclease activity. The nucleases are the main actors of DNA degradation. Some other enzymes, such as neutrophil elastase and myeloperoxidase, that are involved in NETosis, also contribute to cfDNA release22. Earlier research has focused on a limited number of genes through targeted analysis23–26. To date, only three genes DFFB, DNASE1, and DNASE1L3, all of which encode nucleases, have been confirmed to impact cfDNA digestion and fragment characteristics27. Studies on nuclease-deficient mice have revealed their role in cfDNA degradation and shaping fragmentomics11,28,29. It has been demonstrated that cfDNA is first generated intracellularly with DFFB, DNASE1L3, and other nucleases, followed by extracellular fragmentation with circulating DNASE1L3 and DNASE130. Despite these findings, the genome-wide genetic basis and additional key genes involved in regulating cfDNA formation remain unclear. Identifying more genes involved in the cfDNA fragmentation process could significantly enhance our understanding of cfDNA dynamics and their roles in various physiological and pathological conditions. Expanding the genetic landscape of cfDNA fragmentation will provide deeper insights into cfDNA’s functions and its potential applications in disease diagnosis and monitoring.

Genome-wide association study (GWAS) is a research method used to identify genetic loci that are statistically associated with specific traits or diseases31. Over the past two decades, GWAS have successfully identified hundreds of thousands of variants linked to numerous human traits and diseases, including height, weight, cardiovascular diseases, tumors, and autoimmune diseases32–34. As the most comprehensive repository of genotype-trait associations, by May 2024, the GWAS catalog contained nearly 7000 publications and over 600,000 association signals, serving as a valuable resource for exploring the genetic background of phenotypes35.

In previous works, our group has proven that with appropriate algorithms and methodology, cfDNA could be utilized as a resource for high-quality genetic analysis, including variant detection, allele frequency estimation, genetic structure analysis, and GWAS36–41. In this study, we leveraged genetic and phenotypic data from 28,016 Chinese pregnant women to conduct the first genome-wide association study on cfDNA motif frequencies, termed cfGWAS (short for cfDNA-GWAS). We used the same copy of whole-genome sequencing data from cfDNA to extract both genotype information and cell-free DNA features, identifying genetic loci associated with these cfDNA molecular characteristics. A comprehensive overview figure has been provided to illustrate the overall design of this work (Fig. 1).

Key findings from our study include the identification of 15 study-wide significant loci linked to 176 cell-free DNA end motifs, the majority of which were novel associations to cfDNA features, with rediscovery of well-known cfDNA-related genes such as DFFB and DNASE1L3, and the discovery of novel genes like PANX1 and DNASE1L1. These results were validated through independent replication studies and experimental validation in knockout mice and cell lines. Furthermore, our one-sample integrative analysis revealed strong genetic correlations between cfDNA features and phenotypes, such as body mass index (BMI), weight, white blood cell count, and neutrophil count. Given the widespread clinical use of cfDNA whole-genome sequencing, our work underscore the paradigm-shifting potential of cfGWAS in medical research and practice.

Results
Study participants
Participants were recruited from Wuhan Children’s Hospital during routine pregnancy examinations between 2017 and 2020. The only inclusion criterion was the availability of NIPT genotype data. We excluded pregnant women carrying fetuses with trisomies and those with a sequencing depth below 0.05× or a mapping rate under 90%, resulting in a final cohort of 38,668 individuals. All samples were used for genotype imputation to maximize accuracy. Please note that not all 38,668 pregnant women had complete phenotype data, such as age and NIPT test timing. After excluding samples lacking this information, 28,016 samples remained for further analysis. For these 28,016 samples, we provided detailed information, including year of pregnancy, maternal age, gestational week at the time of the NIPT test, type of blood tube used, sex of newborn, and fetal fraction (Supplementary Fig. 1A–F). Approximately 40% of the samples were collected in 2020, with the remainder collected between 2017 and 2019. The average maternal age was 29 years, and the average gestational age at the time of NIPT testing was 16.16 weeks. The male-to-female newborn ratio was 1.10, and the average fetal fraction was 9.8%, consistent with previous reports of 10%42,43.

Genotype and phenotype data
Notably, starting in December 2018, an upgraded NIPT product led to increased sequencing depth. After excluding 1719 samples with unknown collection dates, there were 10,795 pregnant women with an average sequencing depth of 0.1049× before December 2018, and there were 26,154 participants with an average depth of 0.1644× afterward, resulting in an overall average depth of 0.1525× (Supplementary Fig. 1G). The sequencing read length remained consistent at 35 base pairs (bp).

Following STITCH imputation (see the “Methods” section), we obtained a total of 8,134,302 single-nucleotide polymorphisms (SNPs). After filtering out SNPs with an imputation info score <0.4, a Hardy–Weinberg equilibrium (HWE) p-value <1e-6, and a minor allele frequency (MAF) <0.05, 2,897,824 high-quality SNPs remained. To assess imputation accuracy, we calculated Pearson’s correlation between the imputed genotype dosages and the original high-depth genotypes from 30 Chinese samples in the 1000 Genomes Project (1KGP). On average, the imputation accuracy (R) was 90.46%, with an R² of 82.79% (Supplementary Fig. 1H). Additional imputation evaluation results, including analyses using other datasets and SNP-level performance, are detailed in our previous paper36,41. It is notable that the equivalent sample size for GWAS is the product of the true sample size and the imputation accuracy R² value.

Additionally, we calculated the Pearson’s correlation (R²) between the alternate allele frequencies of these high-quality SNPs in our dataset and those in the East Asian samples from the 1000 Genomes Project (Supplementary Fig. 1I). An R² of 0.985 demonstrates a high concordance between our imputed variants and those from the 1KGP.

The phenotypes of interest are cell-free DNA molecular features, specifically the frequencies of 4-mer end motifs. These motifs represent the first four nucleotides at the 5’ end of each plasma DNA strand, leading to a total of 256 distinct phenotypes. We presented the distributions of the 256 motif frequencies using bar plots and box plots (Supplementary Fig. 2). Overall, the distributions of motif frequencies are consistent with previously reported patterns11. For example, the motif CCCA exhibit the highest frequency, C-end motifs tend to have higher average frequencies than others, and motifs containing CG dinucleotides generally show lower frequencies.

GWAS identified 15 highly credible loci
In the GWAS analysis, we included maternal age, gestational week at NIPT screening, sequencing depth, and the first 10 PCs of the genotype data as covariates. Since not all 38,668 participants had complete covariate information, the adjusted sample size was reduced to 28,016. At the genome-wide significance threshold of 5e−08, we identified a total of 507 significant motif-locus associations for the 256 end motif frequencies, involving 32 unique loci and 208 end motifs (Supplementary Data 1). While at the stricter study-wide significance threshold of 1.95e−10 (=5e−8/256), there were 345 significant associations, involving 15 loci and 176 motifs (Supplementary Data 2). We regarded these 15 loci as a group of highly credible signals and conducted follow-up analyses for those not utilizing the complete GWAS summary statistics. To assess the reliability of our GWAS results, we drew the quantile–quantile plots (QQ-plots) and calculated the genomic inflation factor (GIF) of each motif (Supplementary Fig. 3). Specifically, the GIF values were consistently around 1 across all 256 motifs, indicating no substantial inflation and suggesting well-conducted GWAS analyses.

To provide informative visualization of our GWAS results, we synthesized the summary statistics for the 256 motifs by selecting the minimum p-value for each SNP across all motifs, which was used to generate a Manhattan plot (Fig. 2A). The most significant locus identified is PANX1, with the lead SNP rs1138800 (p-value=7.67e−298, beta=3.65e−01, se=9.76e−03). Following PANX1 is DFFB, with lead SNP rs34725762 (p-value=1.82e−286, beta=4.53e−01, se=1.24e−02). The subsequent two most significant loci were DNASE1L3 (rs12633655, p-value=5.13e-30, beta=−1.12e-01, se=9.82e−03) and DNASE1L1 (rs2283762, p-value=8.44e−25, beta=1.95e−01, se=1.89e−02). For more details of each signal, we provided the regional plots (Supplementary Fig. 4). Other loci reaching the study-wide significant threshold include SMIM24, PSMD3, FCHO2, PADI4, ABO, and so on. For each locus, we counted the number of its associated motifs (Fig. 2B). The top five loci with the largest number of associated motifs were PANX1 (136), DFFB (123), PSMD3 (21), DNASE1L3 (17), and FCHO2 (9).

From the perspective of end motifs, we summed up the number of significant loci for motifs with the same starting nucleotide. The counts for ‘A’, ‘T’, ‘C’, and ‘G’ were 90, 97, 80, and 78, respectively.

We then focused on highlighting the biological functions of four key loci: PANX1, DFFB, DNASE1L3, and DNASE1L1. The gene PANX1 encodes Pannexin-1, a cell membrane channel protein that mediates intracellular–extracellular molecule exchange and participates in cell signaling, apoptosis, inflammation, and development44,45. The gene DFFB (DNA fragmentation factor subunit beta) encodes a protein essential for apoptosis. Specifically, DFFB is a key component of the DNA fragmentation complex, responsible for cleaving DNA into fragments during apoptosis30. The gene DNASE1L3 encodes deoxyribonuclease 1-like 3, a nuclease that degrades DNA. It is predominantly expressed in the spleen and lymph nodes, and is crucial for clearing debris from apoptotic cells, thereby preventing autoimmune responses triggered by self-DNA and maintaining tissue homeostasis30. Similarly, the gene DNASE1L1 encodes deoxyribonuclease 1-like 1, another member of the deoxyribonuclease I family. However, no established links exist between DNASE1L1 and diseases caused by deficient chromatin clearance during apoptosis or necrosis in vivo46.

While DFFB and DNASE1L3 are established apoptotic nucleases crucial for cfDNA fragmentation30, we turned to the novel gene PANX1 for further investigation. To this end, we conducted GWAS replication in three independent cohorts and performed experimental validation using both knockout models and knockout cell lines.

Information for GWAS replication studies
Three independent replication cohorts were introduced and utilized to validate the study-wide associated loci identified in this study. In the Chinese NIPT replication cohort, we recruited 25,578 pregnant women who underwent NIPT tests at Wuhan Children’s Hospital between 2020 and 2021. We provided information on the recruitment year and age distribution of these samples (Supplementary Fig. 5A and B). The Dutch replication cohort was obtained from a recent GWAS of NIPT samples47. They performed GWAS analysis for six cfDNA features, including diversity of bincount, fetal concentration, mitochondrial DNA concentration, fragment size diversity, motif diversity score, and total cfDNA concentration. The sample sizes in GWAS analysis for each phenotype are provided in Supplementary Fig. 5C. To validate our GWAS findings in the general population, rather than in pregnant women, we introduced a natural population replication cohort consisting of 442 normal healthy individuals, consisting of 206 males and 261 non-pregnant females, with an average age of 34.8 (Supplementary Fig. 5D and E).

Replication on integrated study-wide associated loci
We first examined the 15 study-wide significant loci identified by integrating GWAS results across all 256 motifs. At a replication significance threshold of 1e−5, 14 (93.3%), 14 (93.3%), and 3 (20.0%) loci were replicated in the Chinese NIPT, Dutch NIPT, and natural population cohorts, respectively (Fig. 2C, Supplementary Data 2). Specifically, PANX1, DFFB, DNASE1L3, and MSR1 were replicated in the Chinese NIPT cohort; PANX1, DFFB, DNASE1L3, and DNASE1L1 in the Dutch NIPT cohort; and PANX1, DFFB, and PADI4 in the natural population cohort. The replication rates were notably higher in the pregnant cohorts than in the non-pregnant cohort, likely due to larger sample sizes. These findings suggest that the identified genetic associations are not limited to the pregnant population but also in the general population.

Overall, three loci (PANX1, DFFB, and PADI4) were consistently replicated across all three studies, indicating shared genetic associations between pregnant and non-pregnant individuals. Additionally, 13 loci (86.7%) were replicated by at least two studies, while all loci were replicated by at least one study. Ten loci—including ABO, MSR1, and FCHO2—were replicated by two NIPT studies but not in the natural population, suggesting they represent pregnancy-specific genetic associations. Notably, FCHO2 is highly expressed in female reproductive organs, including the uterus, vagina, and cervix, as shown in Supplementary Fig. 6 from the GTEx database (https://gtexportal.org/home/).

Replication on individual motif-locus associations
We assessed the replication status of the 345 study-wide significant motif-locus association pairs using the GWAS summary statistics of the corresponding motifs in the Chinese NIPT and natural population replication study. Overall, 256 associations (74.2%) were replicated in the NIPT study and 44 (12.8%) in the natural population study, with a total of 259 unique associations replicated in at least one cohort (Supplementary Data 3, Supplementary Fig. 7). For the 256 and 44 replicated motif-locus associations in the Chinese NIPT and natural population cohorts, respectively, we further examined whether their top SNPs showed consistent effect directions between the discovery and replication studies (Supplementary Fig. 8, Supplementary Data 3). The scatter plots demonstrated that all replicated motif-locus associations (solid dots) exhibited consistent effect directions across the corresponding studies. For each motif, we calculated the Pearson’s correlation coefficients (R) between genetic effects in the discovery study and replication studies (Supplementary Fig. 9). Among the 176 significant motifs, 159 showed R>0.5, and 143 showed R>0.75 in the NIPT cohort, while 132 and 112 met these thresholds in the natural population cohort, respectively. The NIPT cohort exhibited a higher degree of consistency in both the direction and magnitude of genetic effects, whereas the natural population cohort showed lower concordance, likely due to its smaller sample size and the non-pregnant samples.

Of the 345 significant motif-locus associations identified in this study, 341 (98.84%) could be replicated in at least one of the six cfDNA phenotypes in the Dutch cohort (Supplementary Data 3). Notably, the “cleave-site motif diversity” phenotype appeared most frequently in both metrics (Supplementary Fig. 10), suggesting it has the strongest association with cfDNA end motif frequencies among the six features analyzed.

Experimental verification in PANX1 knockout mice and cell lines
To validate the role of PANX1 in cell-free DNA generation and clearance, we conducted knockout (KO) experiments in both mice and cell lines. We compared plasma cfDNA motifs between Panx1 KO mice (n=6) and wild-type (WT) mice (n=6), as well as cfDNA motifs in cell culture supernatants from PANX1 KO and WT Jurkat cell lines (5 replicates each) (Methods).

In mice, the motif diversity score (MDS) of plasma cfDNA was significantly lower in Panx1 KO mice (mean: 0.923, range: 0.916–0.932) compared to WT mice (mean: 0.930, range: 0.928–0.935, t(10)=−2.78, p-value=1.95e−02, mean difference=−0.007, 95% CI=[−0.0128, −0.0014]) (Fig. 3A). This alteration in cfDNA motif profiles after Panx1 knockout supports the gene’s key role in cfDNA fragmentation. We further examined cfDNA motifs associated with PANX1 genetic variants identified in our genome-wide association study. The frequency of motifs positively correlated with PANX1 increased significantly in Panx1 KO mice (mean: 9.71, range: 9.03−10.23) compared to WT mice (mean: 9.11, range: 8.67–9.70, t(10)=2.35, p-value=4.03e−02, mean difference=0.594, 95% CI=[0.0321, 1.1564]) (Fig. 3B). Conversely, the frequency of negatively correlated motifs decreased in Panx1 KO mice (mean: 10.00, range: 8.91–11.31) compared to WT mice (mean: 11.43, range: 10.33–12.35, t(10)=−2.82, p-value=1.83e−02, mean difference=−1.437, 95% CI=[−2.5737, −0.2999]) (Fig. 3C). These findings suggest the impact of PANX1 activity on cfDNA motif patterns.

Similar trends were observed in the PANX1 KO cell line experiments, with an enhanced alteration in motif patterns. The mean MDS decreased from 0.965 in WT cells to 0.956 in KO cells (t(8)=−17.04, p-value=1.43e−07, mean difference=0.009, 95% CI=[0.0078, 0.0102]) (Fig. 3D). The frequency of motifs positively correlated with PANX1 increased significantly in KO cells (mean: 7.02, range: 6.93–7.22) compared to WT cells (mean: 6.10, range: 6.01–6.20, t(8)=13.9, p-value=7.12e−07, mean difference=0.930, 95% CI=[0.7757, 1.0855]) (Fig. 3E). Conversely, the frequency of negatively correlated motifs decreased in KO cells (mean: 16.26, range: 15.91–16.46) compared to WT cells (mean: 18.09, range: 17.83–18.35, t(8)=−14.2, p-value=5.90e−07, mean difference=−1.829, 95% CI=[−2.1261, −1.5319]) (Fig. 3F). These results from both in vivo and in vitro experiments consistently demonstrate the significant role of PANX1 in shaping cfDNA fragmentation patterns.

The fold change of cfDNA concentration in each mouse relative to the mean cfDNA concentration in WT mice was calculated. As depicted in Supplementary Fig. 11, we observe a significant increase in the fold change of cfDNA concentration in Panx1 KO mice (mean: 1.80, range: 1.30–2.90) compared to WT mice (mean: 1.00, range: 0.69–1.43, t(10)=3.08, p-value=1.17e−02, mean difference=0.805, 95% CI=[0.2225, 1.3882]), suggesting the knockout of Panx1 gene influenced the release of cfDNA into plasma (Supplementary Fig. 11A). Although the difference did not reach statistical significance, we do observe a trend of increased cfDNA concentration of PANX1 KO cells (mean fold change: 1.10, range: 1.05–1.15) compared to the WT cells (mean fold change: 1.00, range: 0.86–1.10, t(8)=2.10, p-value=6.89e−02, mean difference=0.098, 95% CI=[−0.0096, 0.2065]) (Supplementary Fig. 11B). The results of both in vivo and in vitro KO experiments suggest that the PANX1 gene plays a crucial role in the generation of cfDNA, and the knockout of the PANX1 gene leads to an increased release of cfDNA.

cfDNA end motif frequencies are heritable
Heritability measures the proportion of phenotypic variance explained by genotypes, thus quantifying the genetic influence on phenotypes. We calculated the heritability of each motif using linkage disequilibrium (LD) score regression (Supplementary Data 4, Supplementary Fig. 12). Notably, among the 256 motifs, 47 motifs exhibited a heritability exceeding 10%, while 78 motifs fell within the range of 5–10%, suggesting a substantial genetic component in their variation.

We also provide the number of genome-wide and study-wide significant loci associated with each motif (Supplementary Data 4, Supplementary Fig. 12). As expected, the heritability is generally positively proportional to the number of associated loci.

Enrichment in hematological cell types
Our heritability partition analysis revealed potentially significant contributions of cell type-specific elements to the heritability of cfDNA end motifs. At a suggestive significance threshold of 1e−2, we identified 323 significant cell-type-motif associations, involving 82 unique cell types and tissues (Supplementary Data 5, Supplementary Table 1, Supplementary Fig. 13A). Interestingly, the most frequently implicated cell types in these associations were parotid gland (33 motifs), followed by salivary glands (17 motifs), neutrophils (17 motifs), blood (16 motifs), and phagocytes (14 motifs), all of which belong to the digestive and blood/immune categories. The most prominent associations are observed between the CGAG motif and two tissues: the parotid gland (p-value=2.57e−5) and the salivary glands (p-value=5.68e−5), and others include the CACC and blood (p-value=1.12e−4), and GTTA and neutrophils (p-value=2.27e−4) (Supplementary Fig. 13B). These findings align with the known role of phagocytes (including neutrophils and monocytes) as both sources and clearers of cfDNA during immune responses48. Our analysis identified an enrichment of hematological cell types, particularly immune cells, in the heritability of cfDNA motif features, indicating their active involvement in cfDNA processing. However, the associations between the parotid gland and salivary gland with cfDNA remain underexplored and lack relevant reports, warranting further in-depth investigation. It should be noted that these results were based on a nominal significance threshold of 0.05 without strict multiple testing correction. Therefore, the findings should be interpreted with caution and regarded as exploratory, providing preliminary indications of potentially associated tissues and cell types.

Pathway-based association analysis
The pathway-based analysis results of all 256 end motifs were provided in Supplementary Data 6. At a suggestive significance threshold of 1e−3, a total of 770 motif-pathway pairs demonstrated significance, and they were distributed among four motif groups starting with ‘A’, ‘T’, ‘C’, and ‘G’, associating with 174, 308, 177, and 111 pathways, respectively (Supplementary Fig. 14A). These significant relationships encompassed 79 unique pathways, categorized into groups such as cell death and apoptosis, immune response, cancer-related pathways, among others (Supplementary Table 2).

Among the identified pathways, the top four with the most significant motifs included the Biocarta SET pathway, leading with 83 associated motifs, followed by the Biocarta Mitochondria pathway (78 motifs), the Biocarta DNA fragment pathway (76 motifs), and Reactome apoptosis-induced DNA fragmentation (50 motifs). Specifically, the Biocarta SET pathway is important in regulating B cells by providing essential survival signals that prevent inappropriate apoptosis. The Biocarta Mitochondria pathway sheds light on how mitochondria contribute to cellular functions like metabolism and apoptosis. Furthermore, both the Biocarta DNA fragment pathway and Reactome apoptosis-induced DNA fragmentation are implicated in DNA breakdown, a pivotal process in apoptosis regulation and cell death.

Moreover, we provided the distribution of pathway-based analysis p-values for all motifs through a box plot, sorted by the median (Supplementary Fig. 14B). Notably, the top pathways are almost consistent with the above observations, including the Biocarta DNA fragment pathway, Biocarta SET pathway, Reactome apoptosis-induced DNA fragmentation, Reactome the NLRP3 inflammasome, and Biocarta Mitochondria pathway. In conclusion, the pathway-based association analysis illuminated the biological processes, such as DNA fragmentation, apoptosis, and the mitochondrial pathway, in governing cfDNA generation and clearance.

Genetic correlation with some pregnancy phenotypes
In a parallel study with the same group of pregnant women, we have performed genetic analysis on their 104 pregnancy phenotypes, including laboratory tests during pregnancy (e.g., hematology, liver function), maternal information (e.g., height, BMI), and neonatal outcomes (e.g., birth weight, birth length)36. We conducted a genetic correlation analysis between the 256 end motifs and the 104 pregnancy phenotypes to assess the overall genetic similarity between each pair (Supplementary Data 7, Fig. 4A). At the strictest Bonferroni-corrected significance threshold of 1.88e−6 (=0.05/(256*104)), we observed a total of 13 significant motif–phenotype pairs, involving maternal BMI and maternal weight. Furthermore, when we relaxed the threshold to a more lenient suggestive level of 1e−3, a total of 138 motif–phenotype pairs showed significant correlations (Supplementary Table 3). To further ensure the reliability of the results, we retained only pregnancy phenotypes with more than five associated motifs, resulting in seven phenotypes (Supplementary Fig. 15): maternal BMI (43 motifs), maternal weight (37 motifs), white blood cell count (13 motifs), neutrophil count (13 motifs), lymphocyte percentage (9 motifs), neutrophil percentage (8 motifs), and uric acid (6 motifs). These genetic correlation results were based on a suggestive significance threshold of 1e−3, not subjected to strict multiple testing correction. Thus, the findings should be interpreted with caution and considered exploratory.

Causal effects of hematological cells on cfDNA motifs
We conducted bi-directional Mendelian randomization analyses to explore the potentially causal relationships between pregnancy phenotypes and end motif frequencies. Initially, we treated pregnancy phenotypes as exposures and motifs as outcomes (denoted as P→M) and subsequently reversed the roles (denoted as M→P). At the Bonferroni-corrected significance threshold of 3.55e−6 based on inverse variance weighted (IVW) method, our P→M analysis identified 168 potentially causal associations involving 10 pregnancy phenotypes, 7 of which belonged to the hematological category; these associations were further confirmed using MRlap, which accounts for sample overlap and other biases (Supplementary Data 8, 9, Fig. 4B). Specifically, 91 out of 168 (54.2%) identified causal relationships were between end motifs and leukocytes, including the white blood cell count (30 causalities), percentage of neutrophils (23), percentage of lymphocytes (19), and neutrophil count (19). In addition, the 57 causal associations between BMI and motifs also constitute a substantial proportion of the results. For motifs starting with ‘A’, ‘T’, ‘C’, and ‘G’, they were causally affected by 36, 67, 48, and 17 phenotypes, respectively.

The relationship between BMI, weight, and cfDNA during pregnancy has been extensively researched. Multiple studies have identified an association between increasing maternal weight and a decrease in cell-free fetal DNA fractions43,49–51. Additionally, other studies have shown that total cfDNA levels are elevated in obese pregnant women, primarily due to increased necrosis and apoptosis of adipose tissue52–54. The relationship between leukocytes and cfDNA has also been extensively studied, which can be elaborated via cellular apoptosis, infection, and inflammatory responses, and DNA clearance mechanisms55,56. The leukocytes, particularly neutrophils, are the major contributors to cfDNA in cancer samples, accounting for ~76% of cfDNA57. Specifically, activated neutrophils undergo NETosis, a specialized form of cell death that forms neutrophil extracellular traps (NETs), significantly increasing cfDNA levels58. Moreover, cfDNA released during exercise mainly originates from extramedullary polymorphonuclear neutrophils, influenced by physical impact, hypoxia, and elevated core body temperature59.

To further ensure the robustness of the 168 causal associations, we performed genetic pleiotropy and heterogeneity tests (Supplementary Data 9). At the significance threshold of 1e−3, all the causal relationships are not affected by genetic pleiotropy, with the smallest p-value being 5.86e−3. We additionally provided the fitted linear graphs for causal pairs of phenotypes and motifs based on the MR-Egger method (Supplementary Fig. 16). The intercepts of the lines were close to zero, also indicating the absence of pleiotropic effects and suggesting reliable causal effects in the MR analysis. Notably, no associations show potential heterogeneity effects (p-value <1e−3) that may impact the overall findings.

We further compared the effect sizes from MR with those from linear regression analysis for the associations between phenotypes and motifs. For both the IVW and MR-Egger methods, the estimated effect directions for the majority of phenotype–motif pairs were consistent with those from regression analysis (Supplementary Fig. 17). The correlation coefficient R of the MR-Egger method is 0.577, and there is also a correlation (R=0.558) between IVW beta and regression beta. This finding highlighted the consistency and differences of these methods.

To ensure the validity of our P→M MR results, we performed reversed M→P analysis, with the cfDNA motifs as exposures and the phenotypes as outcomes. At the Bonferroni-corrected significance threshold of 3.27e−6, we found no causal effects of motifs on phenotypes, with the smallest p-value being 1.28e−5 (TGCC→red cell distribution width-standard deviation). This phenomenon is consistent with our prior knowledge that cfDNA features are often considered as consequences rather than causes; the occurrences of some biological or physiological processes (e.g., immune response) may result in alteration of hematological tests (e.g., neutrophil count), causing the change of cfDNA molecular signatures (e.g., concentration, motifs). Our bi-directional MR analysis, which controlled for both genetic pleiotropy and heterogeneity effects, highlighted the causal effects of hematological indicators, especially immune cells, on cfDNA molecular features, aligning with the understanding that most cfDNA originates from leukocytes.

Colocalized signals between cfDNA motifs and pregnancy phenotypes
To identify shared causal genetic effects on both end motifs and pregnancy phenotypes, we conducted GWAS-GWAS colocalization analysis. This analysis revealed 130 colocalized motif-phenotype pairs at a posterior probability of H4 (one common causal variant) >0.75, distributed across five loci: PSMD3 (94 pairs), FCHO2 (12 pairs), ZFPM2 (9 pairs), ABO (9 pairs), and DEFA10P (6 pairs) (Fig. 4C, Supplementary Fig. 18, Supplementary Data 10). Notably, the phenotypes with the highest number of colocalized signals were primarily related to hematological traits (leukocytes, thrombocytes, and erythrocytes) and protein levels (globulin and albumin).

The most frequent colocalizations appear in locus PSMD3, linking cfDNA motifs and leukocyte cells. The gene PSMD3 (proteasome 26S subunit, ATPase, 3) plays critical roles in the process of protein degradation, in maintaining cellular health, and in regulating various cellular processes. Silencing of PSMD3 inhibits cell proliferation and induces apoptosis, as evidenced by the significant reduction in breast cancer cell proliferation and colony formation following Psmd3 knockdown60. Previous GWAS studies have extensively reported their associations with leukocyte cells61–64 and immune diseases (e.g., asthma, allergic disease)65,66.

Additionally, from the perspective of end motifs, T-end exhibited the largest number of colocalized signals (91), followed by C-end (17), A-end (13), and G-end (9). These findings provide evidence for a shared genetic background influencing both cfDNA end motifs and pregnancy phenotypes. The significant enrichment of hematological and protein-related traits in the colocalized signals suggests their close relationship with cfDNA molecular features.

Pleiotropy analyses with the GWAS catalog
From the GWAS catalog database, we extracted traits that exhibited pleiotropy with the study-wide associated loci reported in this study (Supplementary Data 11). We focused on traits with a high proportion of pleiotropic loci (ratio >4%) as well as those with a large absolute number of pleiotropic loci (count >5) (Supplementary Fig. 19A). In summary, these primarily included a series of hematological traits. Several platelet-related indicators were highlighted due to their high pleiotropy ratios, including platelet forward scatter and platelet counts at different stages of pregnancy. In addition, a range of red and white blood cell-related traits, such as hemoglobin, neutrophil, and lymphocyte, also showed considerable pleiotropy. The genetic overlap between these leukocyte traits and cfDNA motifs provides further genetic evidence for their strong correlation. Furthermore, notable results also included body measurement traits such as height and BMI, as well as total cholesterol levels and low-density lipoprotein cholesterol levels.

Since we did not have cholesterol data for the pregnancy phenotypes, we did not investigate its association with cfDNA motifs. However, previous studies have reported that stressed individuals with cfDNA levels exceeding 1000 showed elevated cholesterol and LDL cholesterol levels67. Additionally, plasma cfDNA has been proposed as a novel biomarker for monitoring and diagnosing cardiovascular diseases68. A recent work based on cfDNA GWAS also reported pleiotropic effects of identified SNPs on cfDNA molecular features and cardiovascular diseases47.

Pleiotropy analyses with OMIM compendium
Among the 50 mapped genes with study-wide significant SNPs associated with cfDNA end motif frequencies, 15 genes across 10 loci have documented disease associations in the OMIM database (Supplementary Table 4, Supplementary Fig. 19B). These associated diseases are distributed across different organ systems, including nervous system and developmental disorders, blood and immune system disorders, and cardiovascular and circulatory system disorders.

Among these diseases, we observed that some are linked to cfDNA characteristics, including systemic lupus erythematosus (SLE, associated with DNASE1L3) and Barrett’s esophagus (BE, associated with MSR1). In detail, SLE is an autoimmune disease characterized by inflammation with the production of autoantibodies. There is increasing evidence that failure in the clearance of cfDNA by deoxyribonucleases (DNASES), particularly DNASE1L3, can lead to the generation of anti-DNA antibodies and SLE69–71. The plasma cfDNA levels are closely related to the disease severity of SLE and can potentially be used as a biomarker for SLE diagnosis. The gene MSR1, macrophage scavenger receptor 1, plays an essential inflammatory role in multiple processes, including immunity, lung and liver disease, and cancer72. Previous studies have reported that mutations in MSR1 are associated with the presence of Barrett’s esophagus73. CfDNA analysis becomes a promising approach for monitoring neoplastic progression in Barrett’s esophagus74. The identification of these associations between motif-associated genes and genetic disorders suggests a complex interplay between cfDNA end motifs, genetic variations, and human health.

Methods
Experimental model and study participant details
Subjects. The participants were recruited from Wuhan Children’s Hospital during their routine pregnancy examinations from 2017 to 2020. A variety of tests were conducted throughout the entire duration of pregnancy, including non-invasive prenatal testing (NIPT), biochemistry assessments, oral glucose tolerance tests, ultrasound screenings, and more. NIPT is a screening test for fetal chromosomal disorders that analyzes cell-free DNA from a blood sample taken from the mother early in pregnancy. The test is typically performed between weeks 12 and 22 of gestation, though the exact timing may vary by country. Before week 12, the fetal fraction is usually too low to detect genetic abnormalities, while after week 22, terminating a pregnancy due to trisomy can pose greater health risks to the mother. Each pregnant participant provided informed consent before enrollment. This study received approval from the Institutional Review Boards (IRB) of both Wuhan Children’s Hospital (2021R062-E03) and the Bioethics and Biosafety of BGI (BGI-IRB 21088-T2). Additionally, authorization was obtained from the National Human Genetic Resources Management Office (Approval No. [2021] CJ2002). Our sole inclusion criterion was the availability of NIPT genotype data. Pregnant women carrying trisomic fetuses, along with those whose sequencing depth was below 0.05 or mapping rate fell below 90%, were excluded, leaving a final cohort of 38,668 individuals. These individuals were then used for genotype imputation.

Pregnancy phenotypes. Throughout the approximately 40-week gestation period, pregnant women undergo a range of routine blood, liver function, kidney function, and more. In our prior investigation, we curated 104 biochemical examination metrics and birth outcome indicators, conducting a comprehensive genome-wide association analysis36. A detailed list providing the information on these phenotypes, including full name, abbreviation, and phenotypic category, is provided in Supplementary Data 7.

Independent cohort for replication study. To validate our GWAS findings on cfDNA features, we conducted replication studies using three independent cohorts: two NIPT cohorts and one natural population. The first replication study involved Chinese NIPT samples collected from Wuhan Children’s Hospital during 2020 and 2021. We included samples without abnormal fetuses, with a sequencing depth >0.05×, a mapping rate exceeding 90%, and explicit information regarding age, sampling year, and gestational week of the NIPT tests. Ultimately, 25,578 pregnant women were retained for GWAS analysis in the replication study. We note that these samples do not overlap with those from the discovery study. The second replication study utilized Dutch NIPT samples47, in which NIPT cfDNA served as the genotype, while cfDNA molecular features were treated as phenotypes. They investigated six phenotypes: bincount diversity, fetal concentration, mitochondrial DNA concentration, fragment size diversity, motif diversity score, and total cfDNA concentration. We downloaded the full summary statistics for these six phenotypes to determine if our significant loci were also identified in their study. The third replication study comprised 442 natural samples, including both males and females, recruited during health examinations at the Dapeng outpatient department in Shenzhen from 2021 to 2022. Each participant provided informed consent before enrollment. This study received ethical approval from the Institutional Review Boards of the Bioethics and Biosafety of BGI (BGI-IRB 21157-T2). The plasma cfDNA was used for whole genome sequencing, using the DNBSEQ platform with paired-end 100 bp mode and an average sequencing depth of ~35×.

Murine models. The animal study was approved by the Institutional Review Boards of Bioethics and Biosafety of BGI (BGI-IRBA24009) and the Institutional Animal Care and Use Committee (IACUC) of Cyagen (GACU24-SY028). The Panx1-deficient mice model (C57BL/6JCya) was created by CRISPR/Cas-mediated genome engineering at Cyagen. The gene knockout mice were first verified by PCR and gel electrophoresis, revealing the target fragment sizes, and then Sanger sequencing of the PCR products was used to confirm the designed deletion of the locus. The sequencing results of the plasma cfDNA of the Panx1-deficient mice were also checked to confirm that the deletion region in the reference genome had no mapping reads (Supplementary Fig. 20). Housing conditions include a 12-h light–dark cycle, controlled temperature and humidity, and pelleted rodent chow. 8-week-old male and female mice were used for the experiment. Sex differences were not analyzed for Panx1-deficient mice.

Cell line models. PANX1 knockout Jurkat cell line model was created by CRISPR/Cas-mediated genome engineering at Cyagen.

Method details
Library preparation and whole-genome sequencing. Peripheral blood (5 mL) was collected from pregnant women using specialized vacuum blood collection tubes. EDTA tubes were transported at low temperatures (2–8°C), with plasma separation completed within 8 h. For ambient temperature circulating tubes (K tubes, G tubes), samples were transported at 6–35°C, and plasma separation was performed within 96 h. The blood was centrifuged at 1600×g for 10 min at 4°C to collect the upper plasma layer, followed by a second centrifugation at 16,000×g for 10 min to further purify the plasma. Free DNA was extracted from the plasma using a nucleic acid extraction kit, and a DNA library was prepared through end repair, adapter ligation, and PCR amplification. After passing quality control, the DNA libraries were pooled and circularized, and DNA nanoballs were generated via rolling circle replication. Single-end sequencing was performed on a BGISEQ-500 sequencer using the combinatorial probe-anchor synthesis method with a read length of 35 bp. Since NIPT testing is typically conducted around 16 weeks of gestation when the fetal DNA content is ~10%, with half of it being identical to maternal DNA, we did not differentiate between maternal and fetal DNA in subsequent analyses. Instead, we treated all DNA as maternal for analysis purposes.

We employed the fastp software79 to conduct quality control analysis on the sequencing data stored in FASTQ format. This involved the removal of adapter sequences and elimination of low-quality sequence fragments (reads). Following this, we utilized BWA80 to align the quality-controlled reads to the hg38 reference genome81, subsequently converting the aligned reads into BAM format and sorting them. Duplicate reads were subsequently eliminated from the sorted BAM files using the Samtools rmdup tool82. Additionally, we utilized the GATK83 Base Recalibrator to recalibrate base quality scores (BQSR) on the sorted BAM files using known site information. The GATK ApplyBQSR tool was then utilized to perform further base quality score recalibration and sorting on the sorted BAM files, thereby generating index files. Finally, Samtools stats was employed to generate comprehensive statistics for the calibrated BAM files.

Variant detection and genotype imputation. Given the ultra-low nature of NIPT sequencing data, traditional analysis workflows designed for high-depth sequencing, such as variant detection, are inadequate. To ensure precise variant discovery, we used BaseVar (https://github.com/ShujiaHuang/basevar), which was developed for calling variants in ultra-low-pass (<1.0×) sequencing data. Additionally, due to the ultra-low depth of the NIPT sequencing data, covering only about 10% of the genome, there was a significant amount of missing genotype information. To ensure the data were suitable for genetic analysis, we applied genotype imputation using STITCH84 (with parameters K=10, nGen=16,000), which provided allele dosages for the imputed genotypes. To assess the accuracy of this imputation, we randomly selected 30 Chinese samples from the high-coverage 1000 Genomes Project (1KGP) dataset as a reference set. These samples were down-sampled to ~0.1× coverage and imputed alongside the original NIPT data. We then calculated the squared Pearson correlation (R²) between the imputed dosages and the high-coverage genotypes to determine imputation accuracy.

CfDNA end motifs. In accordance with the previous definition19, plasma DNA end motifs were identified by analyzing the first four nucleotides (4-mer) at the 5’ ends of plasma DNA fragments aligned to the reference genome. With four nucleotide positions, there are 256 possible 4-mer motifs. We calculated the frequency of each motif and normalized these frequencies by the total number of fragment ends sequenced. To assess the distribution of motif frequencies, we employed a metric called MDS, based on normalized Shannon entropy. MDS is defined by Eq. (1):

MDS = (1/log256) * Σ (Pi * log Pi) (i=1 to 256)

where Pi is the frequency of the ith motif. The MDS value ranges from 0 to 1, where a higher MDS indicates greater diversity and randomness in the motif distribution. In contrast, a lower MDS reflects a more skewed distribution of motif frequencies. Therefore, a decreasing MDS suggests increasing skewness in the frequency distribution of end motifs, while an increasing MDS indicates that motif frequencies are becoming more uniform.

Murine sample collection. Mice were terminally anesthetized by CO2, their eyeballs were removed and blood was collected in EDTA tubes (KANGJIAN). 6 Panx1−/− mice (3 males and 3 females) and 6 WT mice (3 males and 3 females) were used. The blood was first centrifuged at 1600×g for 10 min at 4°C to collect the upper layer of plasma, and then the plasma was centrifuged at 16,000×g for 10 min at 4°C to remove cellular debris. The supernatant plasma samples were stored at −80°C before cfDNA extraction.

Cell line culture. Gene knockout cell lines were confirmed by sequencing. WT and PANX1 knockout cells were cultured in RPMI 1640 medium at 37°C in a humidified atmosphere with 5% CO2. 10% fetal bovine serum was added to the medium. The cell culture supernatant was first centrifuged at 1600×g for 10 min at 4°C, and then centrifuged again at 16,000×g for 10 min at 4°C to remove cells and cellular debris during sample preparation.

Cell-free DNA extraction and quantification. The blood plasma and cell culture supernatants samples were used for cell-free DNA extraction using MagPure Circulating DNA KF Kit (MD5432-02, Magen). 100 μL plasma from each mouse and 500 μL cell culture supernatant for each sample were used for cfDNA extraction. The cfDNA concentrations were measured by Qubit™ 4.0 Fluorometer (Invitrogen) using Qubit dsDNA High Sensitivity Assay Kit (Q32854, Invitrogen).

Statistics and Reproducibility
Genome-wide association study. We utilized the genotype dosage (a quantitative value between 0 and 2) imputed by STITCH as genotype data, and 256 cfDNA end motifs served as phenotypes for the genome-wide association analysis. To control for population stratification, we performed principal component analysis (PCA) on the population genetic variation identified by BaseVar using PLINK285, selecting the top 10 principal components as covariates. Furthermore, maternal age, gestational age at NIPT testing, and sequencing depth were included as additional covariates. Based on the STITCH imputation results, we filtered SNPs with a minor allele frequency (MAF) greater than 0.05, a Hardy–Weinberg equilibrium (HWE) p-value >1e−6, and an imputation info score >0.4. A SNP list was then created, and only the SNPs in this list (using the --extract parameter) were included for GWAS analysis. Additional arguments included --pheno-quantile-normalize and --covar-variance-standardize, which were used to transform motif frequencies to a standard normal distribution and to standardize quantitative covariates to have a mean of zero and variance of one. We set the significance threshold at 5e−8 for genome-wide significance and 1.95e−10 (=5e−8/256) for study-wide significance.

Validation in replication studies. To validate the results of the GWAS analysis, we utilized three independent replication cohorts and conducted validation through two approaches, based on integrated study-wide associated loci and individual motif–locus associations, respectively. We first focused on the integrated loci. The GWAS summary statistics of 256 motif frequencies from the discovery study were merged and integrated by selecting the minimum p-value across all 256 phenotypes. The GWAS summary statistics from the three replication studies underwent the same merging procedure before proceeding with replication validation. Subsequently, we iterated through a total of 345 motif–locus associations across various motifs in the discovery study to determine whether they could be validated in the corresponding motifs in the three replication cohorts.

In the Chinese NIPT cohort, we applied exactly identical SNP filtering and analytical procedures as those used in the discovery study described above, which yielded a GWAS result comprising 2,794,215 SNPs. In the natural population cohort, the cfDNA end motif extraction method was the same as described above. We filtered SNPs with a MAF >0.01, HWE p-value >1e−6, and a genotype missing rate <10% for the GWAS analysis. The first 10 principal components from PCA, along with age and sex, were included as covariates. Gestational age at blood draw and sequencing depth were not included, as these were not NIPT samples. Regarding the Dutch cohort, since the study did not provide GWAS summary statistics on motif frequency but instead included six cfDNA molecular feature phenotypes, each feature was used for replication when verifying each study-wide motif–locus association. Additionally, it was recorded which specific feature the SNPs supporting the validation originated from.

To define a successfully replicated signal, we extended each significantly associated locus identified in the discovery cohort by 500 kb on both sides and examined whether any SNPs within this extended region in the replication cohorts met the following criteria: (1) a p-value <1e−5 and (2) linkage disequilibrium (r² >0.1) with a significant SNP within the corresponding discovery locus. Additionally, in the replication study based on motif frequency GWAS (Chinese NIPT and natural population cohorts), we assessed the consistency of the effect directions for the top SNP in each discovery locus.

Heritability partition. To further understand the genetic architecture of the cfDNA molecular features and the polygenic contributions to the heritability of different genetic components, including cell type-specific elements, we performed a heritability partition analysis. Specifically, to partition heritability from GWAS summary statistics of the 256 end motifs, we used stratified LD score regression (stratified LDSC)86, which can account for linkage disequilibrium. The stratified LDSC converts cell-type/tissue-trait associations into cell-type/tissue-specific SNP functional annotations and estimates the extent to which different cell-type/tissue-specific functional annotations explain trait heritability. The cell type-specific expression data are referred to the Genotype-Tissue Expression (GTEx) project87 and the Frankel ab dataset88. In total, there are 205 tissues and cell types from nine categories, including adipose, blood/immune, cardiovascular, central nervous system, digestive, endocrine, liver, musculoskeletal/connective, and others. The LD reference data is the 1000 Genomes Phase 3 East Asian populations.

Pathway-based analysis. Pathway enrichment analysis assists in highlighting pathways enriched with candidate genes pivotal for biological functions, thereby facilitating the unraveling of mechanistic insights into the underlying traits. Leveraging summary statistics from GWAS conducted on 256 motifs, we performed pathway enrichment analysis using Pathway Scoring Algorithm (PASCAL) software—an intuitive tool designed for gene scoring and pathway analysis from GWAS results89. The databases include Reactome90, Biocarta database (http://www.biocarta.com/), and KEGG91. For each of the 256 motifs, we computed the median pathway score and prioritized pathways in ascending order based on their scores, with special attention to those with higher rankings.

Genetic correlation. Genetic correlation is defined as an informative metric to quantify the overall genetic similarity between two complex traits. In this section, we investigated the genetic correlation between 256 end motif frequencies and 104 pregnancy phenotypes measured in the same group of individuals. Specifically, we used LD score regression to perform the genetic correlation analysis with the reference panel being the East Asian population from the 1000 Genomes Project92.

Mendelian randomization. Mendelian randomization (MR) is a methodological approach used in epidemiology and genetics to investigate causal relationships between exposure factors and outcomes. In this study, we propose performing bi-directional MR analysis using end motifs and pregnancy phenotypes as either exposure or outcome variables. Specifically, we conducted LD clumping to select only independent SNPs within a 10,000 kb window and a clumping r2 threshold of 0.1. LD was calculated based on the East Asian population from the 1000 Genomes Project. To infer causality, we employed the inverse variance weighted (IVW) method93 and the MR-Egger method94 using the R function TwoSampleMR: mr95. It is worth noting that we only selected exposure variables with more than three significant SNPs after LD clumping and harmonization with outcome data to ensure the reliability of the MR analysis.

In the scenario where pregnancy phenotypes are the exposures and motif frequencies are the outcomes (denoted as P->M), we included 55 exposure variables and 256 outcomes, with a significance threshold of 3.55e-6 (=0.05/(55*256)) after Bonferroni correction. In the reverse scenario, where motif frequencies are the exposures and pregnancy phenotypes are the outcomes (denoted as M→P), we included 147 exposures and 104 outcomes, with a significance threshold of 3.27e-6 (=0.05/(147*104)) after Bonferroni correction. We applied the IVW, weighted median, and MR-Egger methods simultaneously, and defined a pair as significant when the IVW results met the significance threshold. Due to sample overlap between the cfDNA motif GWAS and pregnancy phenotype GWAS, and to minimize potential bias in MR-Egger results, we additionally performed the MRlap method96. MRlap is a two-sample MR method that accounts for sample overlap, weak instrument bias, and winner’s curse.

Colocalization analysis. Colocalization analysis is a statistical method used to determine whether two or more traits share the same causal genetic variant(s) within a specific genomic region. This approach is essential when different GWAS identify overlapping signals, suggesting a potential shared genetic architecture between traits. To investigate whether the cfDNA end motifs and pregnancy phenotypes have overlapped genetic effects due to the same causal variant, we performed colocalization analysis using R package coloc97. The prior probability of a SNP is associated with either end motif, or phenotype is set to be 1e−4, and with both traits is 1e−5. For visualizations of the colocalized regions between end motifs and phenotypes, we used R package locuscomparer98.

Pleiotropy analysis with public databases. In this section, we assessed the pleiotropic effects of motif-associated SNPs/genes based on the GWAS catalog35 and OMIM compendium99. First, we downloaded the latest association data from the GWAS catalog (as of September 29, 2025) and performed initial filtering, retaining only significant SNPs at the whole-genome level (p-value <5e−8). We then defined a search region based on each study-wide association locus associated with cfDNA motifs, extending 100 kb upstream and downstream. For each identified search region, we extracted all significant SNPs from the GWAS catalog data and removed duplicates based on the PubMed ID and associated trait to ensure each trait is recorded only once per study. Then, for each GWAS catalog trait, we counted the total number of recorded associated loci and the number of loci shared with end motifs, and calculated the corresponding ratio. Trait categories were matched using the EFO mappings file from the GWAS Catalog. When a trait was mapped to multiple EFO terms, we prioritized categories not labeled as “Other trait”, “Other measurement”, or “NR”.

OMIM (Online Mendelian Inheritance in Man) is a comprehensive database of human genes and genetic disorders. We created a list for all mapped genes of SNPs associated with cfDNA end motifs and searched for their related diseases in the OMIM database.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Zhu and colleagues describe a GWAS of cfDNA end motif frequencies, carried out using NIPT data from ~28k pregnant women. The authors identify a number of novel genetic regions that may be involved in cfDNA biology, and linked these to pregnancy phenotypes using genetic correlation, MR, and co-localisation approaches. Currently the manuscript lacks sufficient detail in terms of description and/or reporting for several of the analyses. I also have concerns about the appropriateness of some of the methods utilised. Despite this, I think the finding of the paper are of interest, and would likely hold if the below points were addressed.\n1. Sequencing depth was very low, with average depth 0.1525X, and many samples with depths lower than this. In this study, the authors used data from 30 Chinese individuals from the 1KGP downsamples to 0.1X to assess imputation quality, and achieved an R2 of 79.34%. In the STITCH paper, the authors found that a higher sequencing depth of 0.3X achieved a correlation between array genotypes and imputed dosages (r2) of around 0.7.\n• What parameters did the authors use for running STITCH, and how were these determined? Specifically, K, the number of ancestral haplotypes, and nGen, the number of generations since the population was founded. It is recommended (https://doi.org/10.1038/ng.3594) to determine these parameters by using an external replication dataset. Were the 30 1KGP individuals used for this purpose? If so, the R2 values may be an overestimate of the imputation quality.\n• Was any filtering done in terms of imputation quality (INFO score threshold), before calculating the accuracy?\n• Imputation quality is affected by the amount of fetal DNA in a sample, and the low number of sequenced reads in NIPT data (https://doi.org/10.1016/j.celrep.2024.114799). The 1KGP individuals will not be impacted by these NIPT-specific factors, so again, the imputation quality may be an overestimate. Are the authors able to assess imputation quality using individuvals with both NIPT data, and higher coverage non-NIPT data? If such data are not available, this should be noted as a limitation.\n2. What were the distributions of the motif scores? Were any trait transformations performed prior to running the GWAS?\n3. In a few instances, the authors have grouped the results by the motif start (A,C,G, or T). Eg the results on page 8, where the number of loci for motifs with the same starting nucleotide are given, or for the MR results in Figure S14. What is the biological rationale for this?\n4. Key Information is missing from the summary statistics for the GWAS. In order to be useful, the effect, and non-effect allele must be explicitly indicated. What are the allele frequencies for the reported SNPs? It may also be helpful to include imputation quality, if no, or a fairly relaxed INFO score filter were applied.\n5. Results for MPCs and MDS phenotypes are not really provided – only summarised by Manhattan plots. Full results for the top SNPs should be included in a supplementary table.\n6. Not enough information is provided in the replication results – for each cohort, there is just a SNP and a p-value. It is not therefore clear exactly what is being “replicated” for each locus. Specfically:\n• In the Chinese NIPT and natural Chinese replication cohorts, are the same motifs significantly associated as those in the discovery sample?\n• For the motif phenotypes, it is not clear whether direction of effect is the same across discovery/replication cohorts.\n• The SNPs reported in the replication cohorts is often different to the sample – presumably these are good proxies (r2 > ??) with the smallest p-value at that locus, but this is not explicitly stated.\n• For the DutchNIPT samples, the phenotypes not motifs, but related phenotypes; again, for each locus please specify which phenotype(s) were associated.\nI think the authors have attempted to show some of these aspects of replication through figure S6; but these plots are too small to be legibile. Instead, I would recommend the authors provide full association results for all motifs at each locus, to allow better interpretation of what truly replicated.\n7. ACAT analysis is not described in the methods. From the results, it appears the authors used ACAT to combine p-values across 256 motifs, for each SNP. Normally, ACAT is used for combining p-values of several SNPs (eg in a gene,region, pathway), for a single trait. I note that ACAT can be a general method for combining p-values, but I’m not sure whether the optimal method in this case? Can the authors provide additional justification/references for using it in this context? Alternatively, a method like C-GWAS may perhaps be more appropriate? https://doi.org/10.1038/s41467-022-35328-9\n8. For all motifs, a series of analyses were carried out using LDSC: heritability estimated; partioned heritability; genetic correlations.\n• For several motifs in Table S3, heritibility is NA. Why is this?\n• For the motifs with non-missing heritaability estimates, the h2 in many cases is very low, and not significantly different to 0 (p<0.05). For the partioned heritability analyses to be sufficiently powered, it is recomended “the dataset analyzed must have a very large sample size and/or large SNP-heritability, and the trait analyzed must be polygenic” (https://doi.org/10.1038/ng.3404). I do not think this applies to a large number of the motifs analysed. Table S4 reports 52,070 p-values from the partioned heritability analyses for all 256 motifs. If the p-values were uniformly distributed under the null hypothesis, we would expect around 521 of these to be p < 1e-2 (ie the p-value threshold used), just by chance. In fact, only 290 p-values meet this threshold. By running the partioned heritability analyses for all motifs, the authors are greatly increasing their multiple testing burden, while obtaining results for several motifs which are not likely to be reliable, due to the low heritabilites. A better approach would be to run these additional anlysis for motifs only where there is evidence of significant SNP-heritability. Furthermore, I’d recommend applying an FDR adjustment (eg Benjamni-Hochberg) to determine significant results, rather than an arbitrirary nominal threshold.\n• The above point also applies to the genetic correlation analyses.\n9. Regarding the Mendelian Randomisation:\n• How were all of the exposures/outcomes selected? Was this based on the genetic correlations, or just the >3 significant SNP requirement?\n• If I understood correctly, the pregnancy phenotype GWASs used for the MR used the same individuals as the motifs GWASs, thereby making this a one-sample MR. MR-Egger is an approach developed for two-sample MR, which in one-sample MR may indicate a significant causal effect more often than it should (https://doi.org/10.1093/ije/dyab084). There is also a OneSampleMR R package which has other estimators, specifically for one sample MR: https://doi.org/10.32614/CRAN.package.OneSampleMR\n• Given the pregnancy phenotypes and motifs were measured in the same indviduals, it would be useful to compare effect estimates from MR with the observed associations between these variables.\n10. The methods overall are generally lacking in some detail:\n• As mentioned in comment 1, more details are needed regarding the imputation. What parameters were used for STITCH? (ie k and nGen ??). Was any filtering done in terms of imputation quality (INFO score threshold?)\n• For the GWAS, were related individuals removed? Association analyses were undertaken uing Plink, which does not account for realtedness, but it is not stated that the cohort was restricted to an unrelated set.\n• Similarly although ancestry PCs were calculated, with 10PCs used as covariates in the GWAS, it is not mentioned whether any exclusions were made based on ancestry?", "reviewer_2": "This article presents the GWAS results for end motifs in a large dataset of pregnant women. The study identifies 23 study-wide signicant loci, several of which are novel. The study is replicated in three independent cohorts. In addition, one of the genes is further validated by knock out experiments in mice and cell lines. This is a very valuable study, providing novel insights in the factors affecting cfDNA fragmentation. Those insights are of utmost importance for the booming field of cfDNA analyses and cfDNA fragmentomics. The article is well written, the figures are clear and systematic.\nMajor comment:\nThe title suggest the genetic basis of cell-free DNA features is studied. However, the study only examines one of the features of the cell free DNA, namely the end motif frequencies and not other aspects e.g. fragment length, jaggedness, ….. Hence, correct the title accordingly to ‘cfGWAS reveal genetic bases of cell-free DNA end motifs’ would be more suitable.\nMinor comments\nThe study is performed using cfDNA from pregnant women, but the data is analysed as if the cfDNA is derived only from the mother. Since about 10% of the cfDNA is derived from the fetus, which has a different genotype, one could consider there to be an effect on the results/analyses. Would this influence the results/power? How would that affect the heritability scores? Can placenta specific signals be misassigned to the mother? It might be valuable to comment/deliberate this aspect?", "reviewer_3": "In this paper, the authors perform GWAS on data from ~28000 NIPT samples. Their main aim is to unravel genome wide associations with the 256 four-nucleotide cfDNA end-motifs. By doing so, they hope to gain more insight into the biology of cfDNA. As they performed single-end sequencing, they could not take into account fragment length association studies. They used one dataset for the primary analysis, and two subsequent non-overlapping datasets to validate their data. These three datasets were from pregnant women. They finally performed an independent study on a small dataset with data from nonpregnant women. As the PANX1 association was considered to be a novel finding, they subsequently performed experimental verification analysis of the biological role of PANX1. Also, the effects of hematological cells on cfDNA end-motifs was studied and a strong causal relationship between neutrophils and end-motifs was noticed.\nThe study is fully based on statistical analysis, but I did not dive into the statistical methods and trust other reviewers with more relevant expertise will do so.\n1. In general, pages 3-20 are main text and the supplemental data, in which the methods are described, consists of an additional 7 pages (text only). With ~400 words per page (estimation), this means the main text of the paper itself is ~7200 words and the supplemental text ~ 2800 words. Moreover, there are 86 references in the main paper and an additional 25 in the supplemental text. I realize the authors present a very extensive study, with a lot of data, but the paper would certainly benefit from less text. The findings of the involvement in cfDNA biology of genes like DNASE1L3 and DFFB is not new, but because of the many different studies discussions on such not-novel findings also take quite some text. Perhaps the authors could focus more on the added value of finding comparable associations with their study and on the novel findings.\n2. The authors mention the identification of the involvement of the PANX1 gene is a novel finding, but Linthorst et al. (their reference 49) previously also identified an association with cfDNA biology of this gene in their GWAS study. Furthermore, the authors obtained the complete summary statistics of the study of Linthorst et al. but as far as I can see the Dutch group is not acknowledged. I suppose the authors have an agreement with the group on this, but if not, it feels unjust.\n3. In the (recent) past many publications appeared dealing with studies on the biology of cfDNA. This is a very interesting topic, but nevertheless, much is still unknown. I am following these studies with much interest, meanwhile thinking about the clinical relevance / consequences. How can these findings be translated to / be of use in daily practice? Apparently, the authors have some ideas on this, as in the last two sentences of the abstract the authors mention “Novel knowledge uncovered by this study promises to revolutionize liquid biopsy technology” and “the potential of this paradigm is enormous”. Also in the last paragraph of the discussion, the authors mention their research “could contribute to the discovery of therapeutic targets for diseases influenced by abnormal DNA concentrations”, and “the significant discoveries will not only revolutionize our understanding of cfDNA biology but also open new avenues for future clinical applications”. I would like to ask the authors to elaborate more on this.\n4. As the authors mention themselves, a limitation of the study is first of all the relatively small initial cohort size, but more importantly the cohorts consisting of data from pregnant women only. Only the data used for the independent study were from non-pregnant individuals, but this was a small cohort size. When replying on comment 3, please take this also into account. When thinking about clinical applications, will these only be applicable to pregnant women? The two other main fields of interest in cfDNA are cancer and transplantation. Will potential “new avenues” also be applicable to these fields, and also to men or only to women? Are differences to be expected? The authors mention that by now worldwide more than 40 million women have undergone NIPT, but what about data from the other fields and data from men?"}, "author_response": "Response: We sincerely appreciate the reviewer's thorough evaluation of our manuscript and the invaluable feedback provided. Below are our responses and corresponding revisions:\n1. Sequencing depth was very low, with average depth 0.1525X, and many samples with depths lower than this. In this study, the authors used data from 30 Chinese individuals from the 1KGP downsamples to 0.1X to assess imputation quality, and achieved an R2 of 79.34%. In the STITCH paper, the authors found that a higher sequencing depth of 0.3X achieved a correlation between array genotypes and imputed dosages (r2) of around 0.7.\n• What parameters did the authors use for running STITCH, and how were these determined? Specifically, K, the number of ancestral haplotypes, and nGen, the number of generations since the population was founded. It is recommended (https://doi.org/10.1038/ng.3594) to determine these parameters by using an external replication dataset. Were the 30 1KGP individuals used for this purpose? If so, the R2 values may be an overestimate of the imputation quality.\nResponse: We thank the reviewer for this comment. We used K = 10 and nGen = 16,000, where nGen was determined based on the formula nGen = 4 × sample size / K, with a sample size of approximately 40,000. To select an appropriate value for K, we conducted an evaluation using around 20,000 individuals (including 30 high-depth WGS samples as true set), testing K values of 5, 8, 10, 15, and 20. The results showed that while imputation accuracy plateaued at K ≥ 10, computational time increased substantially with higher K values. Therefore, we selected K = 10 as a balance between accuracy and efficiency.\nWe have added and highlighted the following content in the Method section:\nTo ensure the data was suitable for genetic analysis, we applied genotype imputation using STITCH (with parameters K=10, nGen=16,000), which provided allele dosages for the imputed genotypes.\n• Was any filtering done in terms of imputation quality (INFO score threshold), before calculating the accuracy?\nResponse: We thank the reviewer for this comment. For the calculation of imputation accuracy and subsequent GWAS analysis, we included only SNPs with a minor allele frequency (MAF > 0.05) and Hardy-Weinberg equilibrium p-value > 1e-6. After applying these filters, all retained SNPs had an INFO score greater than 0.5. The following figure, which illustrates this, was included in our published work (Supplementary Figure S1C, https://www.cell.com/cell-genomics/fulltext/S2666-979X(24)00237-4):\n[editorial note: third party material redacted]\n• Imputation quality is affected by the amount of fetal DNA in a sample, and the low number of sequenced reads in NIPT data (https://doi.org/10.1016/j.celrep.2024.114799). The 1KGP individuals will not be impacted by these NIPT-specific factors, so again, the imputation quality may be an overestimate. Are the authors able to assess imputation quality using individuvals with both NIPT data, and higher coverage non-NIPT data? If such data are not available, this should be noted as a limitation.\nResponse: We thank the reviewer for this comment. Indeed, the unique characteristics of cfDNA, fetal DNA, and short sequencing reads can influence imputation performance, and genomic DNA data from the 1000 Genomes Project (1KGP) may not fully capture these properties. To address this limitation, in our previous work, we collected data from 14 pregnant women with both high-depth cfDNA sequencing (45.39×) and corresponding real NIPT data. The high-depth cfDNA data were used as the ground truth for genotypes. Under this setting, the imputation accuracy reached approximately 85%. We acknowledge that the high proportion of common variants (and higher imputation accuracy) is mainly due to the limited size (14 samples) of the high-depth sequence data, leading to the exclusion of low-frequency variants from the imputation accuracy calculation. The following figure, which illustrates this, was included in our published work (Supplementary Figure S1D, https://www.cell.com/cell-genomics/fulltext/S2666-979X(24)00237-4):\n[editorial note: third party material redacted]\n2. What were the distributions of the motif scores? Were any trait transformations performed prior to running the GWAS?\nResponse: We thank the reviewer for this comment. We provided the distributions of the 256 motif frequencies using bar plots and box plots, as shown in Supplementary Figure S2. To assess whether these distributions follow a normal distribution, we conducted Kolmogorov–Smirnov (K-S) tests. However, due to the high sensitivity of the K-S test in large samples, none of the motif frequencies passed the normality test. As an alternative, we generated Q–Q (quantile–quantile) plots for several representative motifs, including those with the minimum, two randomly-selected medium, and maximum K-S p-values. Q–Q plots of selected 4-mer motif frequencies suggest that most motif frequencies are approximately normally distributed, although mild deviations—particularly in the tails—are observed. Motifs such as ACAA and GCTC show relatively good alignment with the theoretical normal distribution, whereas ACTG and ATTG exhibit heavier tails and deviations, especially in higher quantiles. These patterns reflect the variability in motif frequency distributions across individuals in a large cohort. We have added and highlighted the following content in the Results section (Phenotype data):\nWe presented the distributions of the 256 motif frequencies using bar plots and box plots (Figure S2). Overall, the distributions of motif frequencies are consistent with previously reported patterns. For example, the motif CCCA exhibits the highest frequency, C-end motifs tend to have higher average frequencies than others, and motifs containing CG dinucleotides generally show lower frequencies (https://pmc.ncbi.nlm.nih.gov/articles/PMC10151549/).\nAdditionally, each motif frequency was normalized to a standard normal distribution for GWAS using the --pheno-quantile-normalize option in PLINK2. The corresponding PLINK2 commands used in the analysis are detailed in the Methods section (Genome-wide association study):\nWe filter in SNPs with minor allele frequency (MAF) greater than 0.05 (--maf 0.05), p-value for testing Hardy-Weinberg equilibrium (HWE) greater than 1e-6 (--hwe 1e-6), and genotype missingness less than 10% (--geno 0.1). Additional arguments included --pheno-quantile-normalize and --covar-variance-standardize, which were used to transform motif frequencies to a standard normal distribution and to standardize quantitative covariates to have a mean of zero and variance of one.\n3. In a few instances, the authors have grouped the results by the motif start (A,C,G, or T). Eg the results on page 8, where the number of loci for motifs with the same starting nucleotide are given, or for the MR results in Figure S14. What is the biological rationale for this?\nResponse: We thank the reviewer for this insightful comment. We acknowledge that our manuscript did not sufficiently elaborate on the biological significance of cfDNA end motifs, beyond the brief mention in the Introduction: “End motifs refer to the nucleotide composition at the 5' end of cfDNA, generated during cfDNA digestion. Analyzing these motifs helps trace cfDNA origin and provides insights into physiological and pathological states.” As end motifs are the primary phenotypes in our study, they indeed merit a more thorough introduction. cfDNA end motifs reflect the activity and cleavage preferences of nucleases. For instance, the nucleases DFFB, DNASE1L3, and DNASE1 preferentially generate cfDNA fragments ending in A, C, and T, respectively. Therefore, the abundance of specific end motifs can indicate the activity of corresponding nucleases. In recent years, cfDNA end motifs have emerged as promising biomarkers in liquid biopsy applications, especially in oncology, as their frequencies have been shown to differ between healthy individuals and cancer patients. Moreover, cfDNA end motifs are tissue-specific and can inform on tissue-of-origin, which has important applications in areas such as transplantation monitoring. We have added and highlighted the following content in the Results section (Phenotype data):\nDifferent cfDNA end motifs reflect nuclease activity, pathophysiological processes, and tissue-specific DNA fragmentation patterns. Specifically, DFFB, DNASE1L3, and DNASE1 preferentially generate A-end, C-end, and T-end motifs, respectively. Recent studies have highlighted cfDNA end motifs as emerging biomarkers for oncology and transplantation monitoring due to their frequency shifts under pathological conditions and their tissue-specific profiles.\nAs for our Mendelian randomization analysis, while we observed that T-end motifs had the greatest number of causal associations with pregnancy-related phenotypes, we currently lack direct biological evidence to fully explain this observation. One possibility is that the unique physiological context of pregnancy may influence placenta-specific nuclease activity, hormonal regulation, or maternal immune adaptation, thereby favoring T-end motif production. We believe this finding has potential clinical significance, suggesting that cfDNA end motif profiling could serve as a novel biomarker for pregnancy-related physiological states, immune responses, or complications. We have added and highlighted the following content in the Discussion section:\nCfDNA end motifs have been increasingly explored as biomarkers in cancer detection and transplantation monitoring. Our team have recently uncovered their role in predicting preeclampsia during early pregnancy. Pregnancy also offers a valuable model to study tissue-of-origin dynamics (maternal vs. fetal). To our knowledge, this is the first study to explore the causal relationships between cfDNA end motifs and physiological and biochemical phenotypes in a pregnant population. We found that T-end motifs exhibited the most causal associations and colocalizations with pregnancy phenotypes. This may be related to placenta-specific nuclease activity, hormonal changes, or maternal immune adaptation during pregnancy. While our current analysis cannot fully explain this observation, it highlights a promising area for further investigation into the biological and clinical significance of cfDNA end motif dynamics in pregnancy.\n4. Key Information is missing from the summary statistics for the GWAS. In order to be useful, the effect, and non-effect allele must be explicitly indicated. What are the allele frequencies for the reported SNPs? It may also be helpful to include imputation quality, if no, or a fairly relaxed INFO score filter were applied.\nResponse: We thank the reviewer for this comment. We will release the full summary statistics for all 256 cfDNA end motif frequencies, including key information such as effect and non-effect alleles, allele frequencies, and imputation quality metrics (e.g., INFO scores). We have submitted a data registration application for these GWAS summary statistics to the Ministry of Science and Technology of China. Upon approval, we will upload the full dataset to the GWAS Catalog, where it will be freely available for visualization and download by the research community. For review purposes, we have uploaded the summary statistics for all 256 GWAS results (organized into four files, each containing 64 motifs) to Google Drive (https://drive.google.com/drive/folders/1Um6bPoBgw4Xokl677T0CDjD_8wxUjTvF?usp=sharing).\n5. Results for MPCs and MDS phenotypes are not really provided – only summarised by Manhattan plots. Full results for the top SNPs should be included in a supplementary table.\nResponse: We thank the reviewer for this comment. After careful consideration, we have decided to focus this manuscript on the 256 distinct cfDNA end motifs, which yielded the most statistically significant GWAS associations in our analysis. In contrast, the MPC and MDS phenotypes showed fewer genome-wide significant associations. To maintain clarity and ensure a focused narrative, we have removed these results from both the main text and the supplementary materials. However, we recognize the potential value of these results and will consider presenting them in a future study or data release for researchers interested in further exploration.\n6. Not enough information is provided in the replication results – for each cohort, there is just a SNP and a p-value. It is not therefore clear exactly what is being “replicated” for each locus. Specfically:\n• In the Chinese NIPT and natural Chinese replication cohorts, are the same motifs significantly associated as those in the discovery sample?\nResponse: We thank the reviewer for this comment. To provide a clearer picture of the replication performance, we have added a new supplementary table (Table S3) containing detailed information on the replication results. This includes all 382 study-wide significant motif–locus association pairs identified in the discovery cohort, along with the top associated SNPs in the replication cohorts, their matched motifs, genetic effect sizes, p-values, and linkage disequilibrium with the originally identified SNPs. In addition, we included a bar plot and a Venn diagram (Supplementary Figure S7) to visually summarize the replication results. Among the 382 study-wide significant motif–locus associations identified in the discovery cohort, 241 (63.1%) were replicated in the Chinese NIPT cohort and 46 (12.0%) in the natural Chinese population cohort, with a total of 244 unique associations replicated in at least one cohort. We have added and highlighted the following content in the Results section (GWAS replication studies):\nIn Table S3, we present detailed replication results for the 382 study-wide significant motif–locus association pairs. Using the full GWAS summary statistics of the 256 motif frequencies from both the NIPT and natural population cohorts, we evaluated the replication of each association. For every motif–locus pair, we report the lead SNP identified in the replication cohort, along with its GWAS p-value and linkage disequilibrium information relative to the corresponding discovery SNPs. To visually summarize the findings, we included a bar plot and Venn diagram (Figure S7). Overall, 241 associations (63.1%) were replicated in the NIPT cohort and 46 (12.0%) in the natural population cohort, with a total of 244 unique associations replicated in at least one cohort.\n• For the motif phenotypes, it is not clear whether direction of effect is the same across discovery/replication cohorts.\nResponse: We thank the reviewer for this comment. To assess the consistency of effect directions between the discovery and replication cohorts (NIPT and natural population), we generated scatter plots comparing the genetic effects of significant SNPs associated with the 256 end motif frequencies in the discovery dataset and those of the same SNPs in the replication cohorts. Additionally, we calculated Pearson’s correlation coefficients (R) to quantify the concordance of genetic effects between the studies. Among the 256 motifs, 158 showed R > 0.5 and 146 showed R > 0.75 when comparing the discovery and NIPT cohorts. For the natural population cohort, 69 motifs had R > 0.5 and 23 had R > 0.75. These results indicate a substantial degree of consistency in the direction and magnitude of genetic effects across datasets. We have added and highlighted the following content in the Results section (GWAS replication studies):\nTo further assess the consistency of effect directions between the discovery and these two replication cohorts, we generated scatter plots comparing the genetic effects of significant SNPs associated with each of the 256 end motifs in the discovery GWAS and the corresponding effects in the replication cohorts (Figure S8). For each motif, we also calculated Pearson’s correlation coefficients (R) between genetic effects in the discovery study and replication cohorts (Figure S9). Among the 256 motifs, 158 showed R > 0.5 and 146 showed R > 0.75 in the NIPT cohort, while 69 and 23 met these thresholds in the natural population cohort, respectively. The NIPT cohort exhibited a higher degree of consistency in both the direction and magnitude of genetic effects, whereas the natural population cohort showed lower concordance, likely due to its smaller sample size and the non-pregnant samples.\n• The SNPs reported in the replication cohorts is often different to the sample – presumably these are good proxies (r2 > ??) with the smallest p-value at that locus, but this is not explicitly stated.\nResponse: We thank the reviewer for this comment. To define a successfully replicated signal, we extended each significantly associated locus identified in the discovery cohort by 500 kb on both sides and examined whether any SNPs within this extended region in the replication cohorts met the following criteria: (1) a p-value < 1e-5 and (2) linkage disequilibrium (r² > 0.1) with a significant SNP within the corresponding discovery locus. The r² values between the associated SNPs in the discovery cohort and the lead SNPs in the replication cohorts are provided in Supplementary Table S3. We have added and highlighted these replication processes in the Method section (Independent cohort for replication study):\nTo define a successfully replicated signal, we extended each significantly associated locus identified in the discovery cohort by 500 kb on both sides and examined whether any SNPs within this extended region in the replication cohorts met the following criteria: (1) a p-value < 1e-5 and (2) linkage disequilibrium (r² > 0.1) with a significant SNP within the corresponding discovery locus.\n• For the Dutch NIPT samples, the phenotypes not motifs, but related phenotypes; again, for each locus please specify which phenotype(s) were associated.\nI think the authors have attempted to show some of these aspects of replication through figure S6; but these plots are too small to be legibile. Instead, I would recommend the authors provide full association results for all motifs at each locus, to allow better interpretation of what truly replicated.\nResponse: We thank the reviewer for this valuable comment. In response, we have added a comprehensive supplementary table (Table S3) that provides the full replication results across all cohorts. For the NIPT and natural population replication studies, we included complete association statistics for each of the 382 study-wide significant motif-locus pairs, along with visual summaries such as bar plots, Venn diagrams, and scatter plots.\nRegarding the Dutch cohort, where replication was based on related cfDNA phenotypes rather than end motifs, we specified in the supplementary table which phenotype(s) were associated at each replicated locus (Table S3). For each significant motif-locus pair, we listed all Dutch phenotypes that supported replication, as well as the top associated phenotype (i.e., the one with the smallest p-value). Furthermore, we summarized the involvement of each phenotype using bar plots: one indicating the total number of motif-locus associations each phenotype was involved in, and another indicating how often each phenotype served as the top phenotype (Figure S10). Notably, the “cleave-site motif diversity” phenotype appeared most frequently in both metrics, suggesting it has the strongest associations with cfDNA end motif frequencies among the six Dutch features analyzed.\nThe above paragraph has been incorporated and highlighted in the revised Results section (GWAS replication studies).\n7. ACAT analysis is not described in the methods. From the results, it appears the authors used ACAT to combine p-values across 256 motifs, for each SNP. Normally, ACAT is used for combining p-values of several SNPs (eg in a gene,region, pathway), for a single trait. I note that ACAT can be a general method for combining p-values, but I’m not sure whether the optimal method in this case? Can the authors provide additional justification/references for using it in this context? Alternatively, a method like C-GWAS may perhaps be more appropriate? https://doi.org/10.1038/s41467-022-35328-9\nResponse: We thank the reviewer for this comment. In the Discussion section of the ACTA paper (https://pmc.ncbi.nlm.nih.gov/articles/PMC6407498/), the authors noted that “As an omnibus testing procedure, ACAT in principle can be applied to combine complementary methods in nearly all kinds of genetic studies, including single-variant analysis, multiple-traits analysis, and set-based analysis such as that considered in this paper.” Based on this, we understood that ACAT is applicable to multiple-trait analyses, which aligns with our study involving 256 end motif traits. However, after careful consideration, we have chosen to focus our manuscript on the 256 individual end motifs, which we believe represent the most statistically robust and interpretable GWAS associations in our analysis. Given that the ACAT analysis is not the primary focus of our study, we have removed these results from the main materials to maintain clarity.\n8. For all motifs, a series of analyses were carried out using LDSC: heritability estimated; partioned heritability; genetic correlations.\n• For several motifs in Table S3, heritibility is NA. Why is this?\nResponse: We thank the reviewer for this comment. In supplementary table, the heritability values are marked as \"NA\" for certain motifs because their estimated SNP heritability was negative when using LDSC. While negative heritability is not biologically meaningful, it can arise due to statistical estimation error. Common causes include small sample size, noisy summary statistics, or inadequate control of population stratification. In our study, the sample size and covariate adjustment were consistent across all 256 motifs. Specifically, the same number of samples were used for all motifs, and population structure was controlled for using the top 10 principal components. Since the majority of motifs showed positive heritability estimates, we believe that small sample size and population stratification are unlikely to explain the observed negative estimates. As mentioned in our response to Question 2, we conducted Kolmogorov–Smirnov (K–S) tests to assess the normality of motif frequency distributions. To explore whether deviations from normality might be related to negative heritability estimates, we generated a scatter plot of K–S test statistics (used instead of p-values due to their near-zero values) against the raw heritability values. However, we did not observe a significant correlation between the two (p = 0.57, linear regression). At this stage, we acknowledge that we have not identified a definitive reason why a small subset of motifs yielded negative heritability estimates. We have added a statement discussing this observation in the Limitations section of the manuscript.\nMoreover, we observed that a small subset of motifs yielded negative SNP heritability estimates using LD Score Regression. While this is likely due to statistical noise rather than biological factors, the exact cause remains unclear and warrants further investigation in larger cohorts or with alternative methods.\n• For the motifs with non-missing heritaability estimates, the h2 in many cases is very low, and not significantly different to 0 (p<0.05). For the partioned heritability analyses to be sufficiently powered, it is recomended “the dataset analyzed must have a very large sample size and/or large SNP-heritability, and the trait analyzed must be polygenic” (https://doi.org/10.1038/ng.3404). I do not think this applies to a large number of the motifs analysed. Table S4 reports 52,070 p-values from the partioned heritability analyses for all 256 motifs. If the p-values were uniformly distributed under the null hypothesis, we would expect around 521 of these to be p < 1e-2 (ie the p-value threshold used), just by chance. In fact, only 290 p-values meet this threshold. By running the partioned heritability analyses for all motifs, the authors are greatly increasing their multiple testing burden, while obtaining results for several motifs which are not likely to be reliable, due to the low heritabilites. A better approach would be to run these additional anlysis for motifs only where there is evidence of significant SNP-heritability. Furthermore, I’d recommend applying an FDR adjustment (eg Benjamni-Hochberg) to determine significant results, rather than an arbitrirary nominal threshold.\nResponse: We thank the reviewer for this insightful comment. We fully agree that the relatively low SNP heritability estimates for many motifs limit the statistical power of the partitioned heritability analyses, as highlighted in the referenced study (https://doi.org/10.1038/ng.3404). As the reviewer pointed out, if the p-values were uniformly distributed under the null hypothesis, we would expect approximately 521 p-values < 0.01, whereas only 290 such p-values were observed, suggesting limited signal. We had indeed attempted several multiple testing correction approaches, including Bonferroni correction and Benjamini–Hochberg (BH) FDR correction. However, no associations remained significant after correction. Despite this, we believed that reporting nominally significant results (p < 0.05) could still provide useful hints about potential motif–cell type associations, while clearly acknowledging their exploratory nature. In light of the reviewer's suggestion, we have explicitly stated this limitation in the Results section (\"Enrichment in hematological cell types\"):\nIt should be noted that these results were based on a nominal significance threshold of 0.05 without strict multiple testing correction. Therefore, the findings should be interpreted with caution and regarded as exploratory, providing preliminary indications of potentially associated tissues and cell types.\n• The above point also applies to the genetic correlation analyses.\nResponse: We thank the reviewer for this insightful comment. As we stated in the manuscript, \"At the strictest Bonferroni-corrected significance threshold of 1.88e-6 (=0.05/ (256*104)), no motif-phenotype pairs reached significance. However, when we relaxed the threshold to a more lenient suggestive level of 1e-3, a total of 61 motif-phenotype pairs showed significant correlations (Table S9). To further ensure the reliability of the results, we retained only pregnancy phenotypes with more than five associated motifs.\" As the reviewer pointed out, the low SNP-heritability estimates for many motifs also limit the statistical power of the genetic correlation analyses, resulting in no motif-phenotype pairs passing the stringent multiple testing correction threshold. Additionally, we would like to clarify that LDSC produces p-values with a minimum resolution of 1e-4; therefore, for four pairs with reported p-values of 1e-4, the true p-values might be even smaller, but the precision of LDSC is limited at this level. We have now added and highlighted the following limitation at the end of this section:\nThese genetic correlation results were based on a suggestive significance threshold of 1e-3 without strict multiple testing correction. Thus, the findings should be interpreted with caution and considered exploratory.\n9. Regarding the Mendelian Randomisation:\n• How were all of the exposures/outcomes selected? Was this based on the genetic correlations, or just the >3 significant SNP requirement?\nResponse: We thank the reviewer for this comment. As noted in the previous response, although some phenotype–motif pairs showed high genetic correlation coefficients (rg), none passed the study-wide significance threshold. Therefore, the MR analyses were not based on the genetic correlation results. Instead, we systematically tested the causality of all phenotype–motif pairs. For the selection of exposure phenotypes, we included only those with more than three independent significant SNPs after LD clumping, ensuring adequate instrument strength. This criterion was applied consistently in both directions of the MR analysis (i.e., phenotype-to-motif and motif-to-phenotype). For the outcome variables, no additional selection criteria were applied—as long as the selected instrumental variables were available in the outcome GWAS dataset and not significantly associated with the outcome, the outcome phenotype was retained for analysis.\n• If I understood correctly, the pregnancy phenotype GWASs used for the MR used the same individuals as the motifs GWASs, thereby making this a one-sample MR. MR-Egger is an approach developed for two-sample MR, which in one-sample MR may indicate a significant causal effect more often than it should (https://doi.org/10.1093/ije/dyab084). There is also a OneSampleMR R package which has other estimators, specifically for one sample MR: https://doi.org/10.32614/CRAN.package.OneSampleMR\nResponse: We thank the reviewer for this insightful comment. We fully agree that the MR-Egger method, while useful for detecting and adjusting for horizontal pleiotropy, can yield biased estimates in the context of one-sample MR due to sample overlap, as highlighted in the cited study (https://doi.org/10.1093/ije/dyab084). To mitigate this, we complemented MR-Egger with the IVW method and only considered results as robustly causal when both methods reached study-wide significance (p < 5e-6). While MR-Egger is sensitive to bias in one-sample settings, its results remain informative if the variability in instrument strength is high. Meanwhile, the IVW method is generally robust in overlapping-sample scenarios and provides a more conservative estimate.\nWe also appreciate the reviewer’s suggestion regarding the OneSampleMR package. Although the exposure (pregnancy phenotypes) and outcome (cfDNA motif frequencies) datasets partially overlap rather than representing a true one-sample setting, we agree that additional validation is necessary. Therefore, we applied the MRlap method (https://github.com/n-mounier/MRlap), which is specifically designed for two-sample MR with potentially overlapping samples and accounts for multiple sources of bias, including weak instruments, winner’s curse, and sample overlap. We re-analyzed the 156 associations that were significant in both MR-Egger and IVW using MRlap. The MRlap results remained consistent, yielding significant associations at the same Bonferroni-corrected threshold (p < 5e-6), with p-value magnitudes comparable to those obtained from the IVW method. These results have been added to Supplementary Table S12. We have included and highlighted the following statements in the Methods section (Mendelian randomization):\nDue to sample overlap between the cfDNA motif GWAS and pregnancy phenotype GWAS, and to minimize potential bias in MR-Egger results, we additionally performed the MRlap method. MRlap is a two-sample MR method that accounts for sample overlap, weak instrument bias, and winner’s curse. We applied MRlap to the 156 associations identified by both MR-Egger and IVW, and confirmed their significance.\nAnd in the Results section (Post-GWAS analyses with one-sample pregnancy phenotypes):\nAt the Bonferroni-corrected significance threshold of 5e-6 based on both MR-Egger method and inverse variance weighted method, our P->M analysis identified 156 potentially causal associations involving 5 pregnancy phenotypes, all belonging to the hematological category; these associations were further confirmed using MRlap, which accounts for sample overlap and other biases.\n• Given the pregnancy phenotypes and motifs were measured in the same indviduals, it would be useful to compare effect estimates from MR with the observed associations between these variables.\nResponse: We thank the reviewer for this insightful comment. In response, we performed linear regression analyses regressing motif frequencies on pregnancy phenotypes for the 156 motif-phenotype pairs identified as significant in the MR analyses. We found that 141 pairs (90.38%) exhibited concordant directions of effect estimates between the MR results and the observed phenotypic associations. Furthermore, the Pearson’s correlation coefficients (R) between the MR-predicted and observed effect estimates were 0.747 for the IVW method and 0.755 for the MR-Egger method, indicating a high degree of consistency.\n10. The methods overall are generally lacking in some detail:\n• As mentioned in comment 1, more details are needed regarding the imputation. What parameters were used for STITCH? (ie k and nGen ??). Was any filtering done in terms of imputation quality (INFO score threshold?)\nResponse: We thank the reviewer for this comment. As detailed in our response to Question 1, we have provided a comprehensive explanation for the selection of STITCH parameters, including the values used for K and nGen. These values have also been added to the revised Methods section for clarity. Regarding imputation quality control, we did not apply an explicit INFO score threshold because, after filtering SNPs with MAF > 0.05 and Hardy–Weinberg equilibrium p-value > 1e-6, all retained SNPs had INFO scores greater than 0.5—commonly considered an acceptable benchmark for downstream analyses.\n• For the GWAS, were related individuals removed? Association analyses were undertaken uing Plink, which does not account for realtedness, but it is not stated that the cohort was restricted to an unrelated set.\nResponse: We thank the reviewer for this valuable comment. In our study, we did not remove related individuals from the cohort. We generated two types of variant files at the population level: one was based on variants called by BaseVar, and the other on variants imputed by STITCH. The BaseVar-called variant files exhibited extensive missingness and were predominantly haplotypic in nature; using these files, the KING software outputted NA values, and PLINK produced -INF kinship estimates, making relatedness inference infeasible. We also attempted kinship estimation based on the STITCH-imputed variants. However, both KING and PLINK reported an unexpectedly large number of related individuals, which we considered likely to be artifacts, given that our cohort primarily consisted of unrelated pregnant women aged between 20 and 40 years.\nTo assess the potential impact of relatedness on our GWAS results, we randomly selected 10 motif frequency phenotypes (covering a range of signal strengths) and reran the GWAS after removing identical twins and first-degree relatives identified above. The p-values were found to be highly consistent between analyses with and without relatedness filtering. This suggests that, given the relatively small proportion of potentially related individuals compared to the overall sample size (>20,000), the presence of relatedness had minimal influence on the results.\nWe have added and highlighted the following content in the Limitation section:\nCurrently, no robust approaches are available for inferring kinship relationships based on NIPT-like ultra-low-depth sequencing data. Neither BaseVar-called variants nor STITCH-imputed genotypes yielded reliable kinship estimates. Given that the cohort consisted predominantly of middle-aged pregnant women, we expect the proportion of close relatives to be very small, and their influence on the GWAS results to be negligible. Nevertheless, the development of appropriate methods for kinship estimation under ultra-low sequencing depth remains an important area for future research.\n• Similarly although ancestry PCs were calculated, with 10PCs used as covariates in the GWAS, it is not mentioned whether any exclusions were made based on ancestry?\nResponse: We thank the reviewer for this comment. Based on the PCA results, we assessed ancestry outliers by examining whether each sample fell within six standard deviations from the mean along both PC1 and PC2. We found no outliers under this criterion and therefore did not exclude any samples based on ancestry. In our previous work (Figure S3d, https://pmc.ncbi.nlm.nih.gov/articles/PMC11602577/), we presented the PC1–PC2 plot, which showed that all NIPT samples clustered centrally, consistent with the recruitment site in Wuhan, China.\n[editorial note: third party material redacted]\nReviewer #2 (Remarks to the Author):\nThis article presents the GWAS results for end motifs in a large dataset of pregnant women. The study identifies 23 study-wide signicant loci, several of which are novel. The study is replicated in three independent cohorts. In addition, one of the genes is further validated by knock out experiments in mice and cell lines. This is a very valuable study, providing novel insights in the factors affecting cfDNA fragmentation. Those insights are of utmost importance for the booming field of cfDNA analyses and cfDNA fragmentomics. The article is well written, the figures are clear and systematic.\nResponse: We sincerely thank the reviewer for the positive and encouraging comments on our work. We greatly appreciate the time and effort the reviewer dedicated to carefully evaluating our study and providing thoughtful feedback. Below are our responses and corresponding revisions:\nMajor comment:\nThe title suggest the genetic basis of cell-free DNA features is studied. However, the study only examines one of the features of the cell free DNA, namely the end motif frequencies and not other aspects e.g. fragment length, jaggedness, ….. Hence, correct the title accordingly to ‘cfGWAS reveal genetic bases of cell-free DNA end motifs’ would be more suitable.\nResponse: We thank the reviewer for this comment. Following the reviewer’s suggestion, we have revised the title to “cfGWAS reveal genetic basis of cell-free DNA end motifs” to more accurately reflect the scope of our study.\nMinor comments\nThe study is performed using cfDNA from pregnant women, but the data is analysed as if the cfDNA is derived only from the mother. Since about 10% of the cfDNA is derived from the fetus, which has a different genotype, one could consider there to be an effect on the results/analyses. Would this influence the results/power? How would that affect the heritability scores? Can placenta specific signals be misassigned to the mother? It might be valuable to comment/deliberate this aspect?\nResponse: We thank the reviewer for this important comment. Indeed, during the NIPT testing window (gestational weeks 12–22), the fetal fraction typically ranges from 5% to 10%. Considering that approximately half of the fetal DNA is maternally inherited, the effective proportion of \"non-maternal\" DNA contamination is estimated to be around 2.5%–5%. In our previous work, we demonstrated that jointly modeling maternal and fetal genotypes during imputation is theoretically feasible and can improve the accuracy of maternal genotype imputation while simultaneously estimating fetal genotype dosages (https://pmc.ncbi.nlm.nih.gov/articles/PMC11602596/). However, we acknowledge that treating all cfDNA as maternal in downstream statistical analyses, without distinguishing the fetal contribution, could introduce some degree of bias. We have now added and highlighted this point as a limitation of the study, as follows:\nNotably, at the time of NIPT testing, the fetal DNA fraction in maternal plasma is approximately 10%, with half of it being identical to the maternal DNA. In our genotype imputation and GWAS analyses, we did not distinguish between maternal and fetal DNA, but rather treated all cfDNA as maternal. While this approach simplifies the analysis, it may introduce minor biases, particularly for loci with strong placental or fetal-specific signals. Recently, a newly developed genotype imputation method, QUITL2, which allows simultaneous imputation of maternal and fetal genomes. In ongoing work, we are applying this method to our NIPT dataset to separate maternal and fetal genomes for future GWAS and post-GWAS analyses, which will help further refine the findings.\nReviewer #3 (Remarks to the Author):\nIn this paper, the authors perform GWAS on data from ~28000 NIPT samples. Their main aim is to unravel genome wide associations with the 256 four-nucleotide cfDNA end-motifs. By doing so, they hope to gain more insight into the biology of cfDNA. As they performed single-end sequencing, they could not take into account fragment length association studies. They used one dataset for the primary analysis, and two subsequent non-overlapping datasets to validate their data. These three datasets were from pregnant women. They finally performed an independent study on a small dataset with data from nonpregnant women. As the PANX1 association was considered to be a novel finding, they subsequently performed experimental verification analysis of the biological role of PANX1. Also, the effects of hematological cells on cfDNA end-motifs was studied and a strong causal relationship between neutrophils and end-motifs and neutrophils was noticed.\nThe study is fully based on statistical analysis, but I did not dive into the statistical methods and trust other reviewers with more relevant expertise will do so.\nResponse: We sincerely thank the reviewer for taking the time to carefully review our manuscript and provide invaluable feedback. Below are our responses and corresponding revisions:\n1. In general, pages 3-20 are main text and the supplemental data, in which the methods are described, consists of an additional 7 pages (text only). With ~400 words per page (estimation), this means the main text of the paper itself is ~7200 words and the supplemental text ~ 2800 words. Moreover, there are 86 references in the main paper and an additional 25 in the supplemental text. I realize the authors present a very extensive study, with a lot of data, but the paper would certainly benefit from less text. The findings of the involvement in cfDNA biology of genes like DNASE1L3 and DFFB is not new, but because of the many different studies discussions on such not-novel findings also take quite some text. Perhaps the authors could focus more on the added value of finding comparable associations with their study and on the novel findings.\nResponse: We thank the reviewer for this valuable comment. After careful consideration, we have decided to focus our manuscript on the analysis of the 256 distinct end motifs, which we believe represent the most statistically robust GWAS associations. Accordingly, we have removed the sections related to MPCs, MDS, and ACAT analyses from both the main text and supplementary materials to enhance the clarity and focus of the manuscript. In addition, we have condensed and streamlined the descriptions throughout the manuscript to improve conciseness. However, we would like to note that due to the incorporation of additional analyses and expanded discussions as requested by the reviewers, the overall length of the manuscript remains relatively substantial.\n2. The authors mention the identification of the involvement of the PANX1 gene is a novel finding, but Linthorst et al. (their reference 49) previously also identified an association with cfDNA biology of this gene in their GWAS study. Furthermore, the authors obtained the complete summary statistics of the study of Linthorst et al. but as far as I can see the Dutch group is not acknowledged. I suppose the authors have an agreement with the group on this, but if not, it feels unjust.\nResponse: We thank the reviewer for this valuable comment. In the study by Linthorst et al., six cfDNA features were analyzed, including bincount diversity, fetal concentration, mitochondrial DNA concentration, fragment size diversity, motif diversity score, and total cfDNA concentration; however, cfDNA end motif frequencies were not among the phenotypes they investigated. We have carefully reviewed our statement and clarified that the association between the PANX1 gene and cfDNA end motif frequencies was newly identified in our analysis. Meanwhile, we also examined the correlation between the six previously studied cfDNA features and the cfDNA end motif frequencies, suggesting that if Linthorst et al. had also analyzed end motif frequencies, PANX1 might have been implicated in their study as well. We have properly cited the work of Linthorst et al. in our manuscript, and their summary statistics are publicly available. In addition, our study represents the first cfDNA end motif GWAS conducted in an Asian population, providing a valuable complement to the Dutch study in terms of both phenotype focus and population diversity.\n3. In the (recent) past many publications appeared dealing with studies on the biology of cfDNA. This is a very interesting topic, but nevertheless, much is still unknown. I am following these studies with much interest, meanwhile thinking about the clinical relevance / consequences. How can these findings be translated to / be of use in daily practice? Apparently, the authors have some ideas on this, as in the last two sentences of the abstract the authors mention “Novel knowledge uncovered by this study promises to revolutionize liquid biopsy technology” and “the potential of this paradigm is enormous”. Also in the last paragraph of the discussion, the authors mention their research “could contribute to the discovery of therapeutic targets for diseases influenced by abnormal DNA concentrations”, and “the significant discoveries will not only revolutionize our understanding of cfDNA biology but also open new avenues for future clinical applications”. I would like to ask the authors to elaborate more on this.\nResponse: We thank the reviewer for this thoughtful comment. Our study primarily focuses on statistical genetic analyses to identify genetic variants associated with cfDNA end motif frequencies, as well as the phenotypic correlations and causal relationships between cfDNA end motifs and pregnancy-related traits. Although we performed experimental validations in knockout mice and cell lines, confirming that cfDNA motif frequencies differ between wild-type and knockout conditions, we did not fully elucidate the underlying biological pathways or mechanisms leading to these differences. Nevertheless, we believe that our findings provide important genetic insights into cfDNA generation and clearance processes. As we discussed in the manuscript, with large sample sizes, cfGWAS has the potential to identify more credible genes associated with cfDNA molecular features, thereby enhancing our understanding of cfDNA biology. This knowledge could ultimately contribute to future clinical applications of cfDNA, such as improving liquid biopsy technologies. For instance, since circulating tumor DNA (ctDNA) often exists in low concentrations, strategies aimed at protecting cfDNA from degradation or reducing its clearance could help maintain sufficient ctDNA levels, allowing sensitive detection even with minimal blood volumes. In addition, we have expanded our discussion regarding the biological significance of cfDNA end motifs. In the revised manuscript, we included the following content:\nIn the Result section:\nWe presented the distributions of the 256 motif frequencies using bar plots and box plots (Figure S2). Overall, the distributions of motif frequencies are consistent with previously reported patterns. For example, the motif CCCA exhibits the highest frequency, C-end motifs tend to have higher average frequencies than others, and motifs containing CG dinucleotides generally show lower frequencies. Different cfDNA end motifs reflect nuclease activity, pathophysiological processes, and tissue-specific DNA fragmentation patterns. Specifically, DFFB, DNASE1L3, and DNASE1 preferentially generate A-end, C-end, and T-end motifs, respectively. Recent studies have highlighted cfDNA end motifs as emerging biomarkers for oncology and transplantation monitoring due to their frequency shifts under pathological conditions and their tissue-specific profiles.\nIn the Discussion section:\nCfDNA end motifs have been increasingly explored as biomarkers in cancer detection and transplantation monitoring. Our team have recently uncovered their role in predicting preeclampsia during early pregnancy. Pregnancy also offers a valuable model to study tissue-of-origin dynamics (maternal vs. fetal). To our knowledge, this is the first study to explore the causal relationships between cfDNA end motifs and physiological and biochemical phenotypes in a pregnant population. We found that T-end motifs exhibited the most causal associations and colocalizations with pregnancy phenotypes. This may be related to placenta-specific nuclease activity, hormonal changes, or maternal immune adaptation during pregnancy. While our current analysis cannot fully explain this observation, it highlights a promising area for further investigation into the biological and clinical significance of cfDNA end motif dynamics in pregnancy.\n4. As the authors mention themselves, a limitation of the study is first of all the relatively small initial cohort size, but more importantly the cohorts consisting of data from pregnant women only. Only the data used for the independent study were from non-pregnant individuals, but this was a small cohort size. When replying on comment 3, please take this also into account. When thinking about clinical applications, will these only be applicable to pregnant women? The two other main fields of interest in cfDNA are cancer and transplantation. Will potential “new avenues” also be applicable to these fields, and also to men or only to women? Are differences to be expected? The authors mention that by now worldwide more than 40 million women have undergone NIPT, but what about data from the other fields and data from men?\nResponse: We thank the reviewer for raising this important point. In our study, we introduced an independent validation cohort composed of non-pregnant individuals and successfully replicated several significant signals such as PANX1, ABO, and MSR1, suggesting that there may be a shared genetic background influencing cfDNA end motif frequencies across both pregnant and non-pregnant populations. At the same time, we also identified signals that were replicated only in the pregnant cohort but not in the natural population cohort, including FCHO2 and PSMD3. Further analysis showed that FCHO2 is highly expressed in female reproductive organs such as the uterus, vagina, and cervix, indicating that it may represent a pregnancy-specific cfDNA signature. As described in the manuscript:\n\"Four loci—FCHO2, PSMD3, and two long non-coding RNA genes—were replicated by two NIPT studies but not in the natural population, whereas NAMPTP2 was replicated only in the natural population and not in the two NIPT studies. Based on these results, we can conclude that the majority of cfDNA-associated genes are shared between pregnant and non-pregnant individuals. If we define pregnancy-specific loci as those replicated by NIPT studies but not in the natural population, then there are four such loci. Notably, FCHO2 is highly expressed in female reproductive organs, including the uterus, vagina, and cervix.\"\nWe acknowledge, however, that the sample size of our non-pregnant validation cohort (~400 individuals) was relatively small compared to the tens of thousands of cases in the pregnant discovery cohorts, leading to lower statistical power in the non-pregnant group. We have clearly stated this limitation in the manuscript: \"Furthermore, our discovery set relied on a cohort of pregnant women, so it’s unclear how many findings are general for all people or specific to pregnancy. Large-scale independent studies on non-pregnant cohorts are necessary to resolve this uncertainty.\"\nRegarding clinical applications, while some signals may be pregnancy-specific, the majority of genetic associations we uncovered are likely to be relevant across different populations, including men and non-pregnant women. Therefore, we believe that the “new avenues” mentioned in our discussion could potentially extend to other major cfDNA research fields, such as oncology and transplantation medicine, although further studies will be necessary to comprehensively characterize any differences."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors provided summary statistics with their revised manuscript. The information provided in these files raises serious concerns regarding the results of the GWAS. While the authors state \"all retained SNPs had INFO scores greater than 0.5\", in these summary results files, there are many SNPs where the INFO_SCORE column has a value < 0.5. A quick look-up in the uploaded file cfGWAS.m256.T___.csv of a couple of the reported SNPs in Table S1 found at least one with very low imputation quality (e.g., rs529723778 has INFO_SCORE of 0.027). I haven't done a comprehensive look-up, but would assume this may be the case for others. That isn't to say that all reported variants are likely to be wrong, but there may be a number of false positive hits reported, and the genome-wide results used in their genetic correlations etc. are likely to contain a lot of noise.\nUnless the INFO_SCORE column in the uploaded summary stats is incorrect, or I have misinterpreted (hopefully this is the case!), much of the analyses should be repeated after filtering out variants with very low imputation quality.\nI address the reviewers' responses to my original comments below:\n1. Overall, the additional detail regarding imputation has adequately addressed my comments. In the response about INFO score thresholds, the authors provide a plot (S1C of their previous paper). It is unclear what the plotted INFO_SCORE is showing here; I'd assume median/mean(?) INFO score, per MAF bin. So while the average INFO score across all SNPs was >0.5, it isn't clear that each SNP had INFO > 0.5. Indeed, the summary statistics uploaded by the authors seemed to confirm that many variants had INFO < 0.5.\n2. Thank you for the additional clarification regarding trait transformation.\n3. The additional explanation is appreciated.\n4. The extra fields requested are available now in the uploaded summary stats but should also be included in the paper (e.g., table S1).\n5. ..\n6. The additional information is welcome. Betas should be provided in table S3 in addition to p-values so consistency of direction of effect can be seen. It appears from Figures S8 and S9 that some variants have opposite directions of effect in the replication cohorts, compared to the discovery. If directions of effects are inconsistent, these shouldn't really be classified as a \"replicated\" association.\n7. ..\n8. ..\n9. The comparison of MR vs observed observations is interesting. Thank you for providing this; perhaps worth considering for a supplementary figure!\n10. Perhaps worth trying to repeat relatedness estimation after filtering for poorly imputed variants.", "reviewer_2": "I am happy with the answers to my comments.\nJust a minor comment:\nIn you answer you state 'QUITL2, which allows simultaneous imputation of maternal and fetal genomes'.\nI believe this method is QUILT2.", "reviewer_3": "Even though the paper was rejected in the first place, the authors submitted a revised version. The authors did a tremendous piece of work, and thoroughly revised the paper. I would like to thank the authors for their efforts!\nAs far as my previous comments are concerned: the authors responded sufficiently to these. My first comment, however, is, despite the response of the authors, not completely resolved, but I leave it to the editor how to proceed with this. The paper is well-written, but it is still very lengthy. Moreover, the results section is full of background information and partly even discussion. I understand why the authors have chosen this format, as otherwise the rationale behind the different analyses would not be understood and the reader would be completely lost, but the fact is that this makes the paper very lengthy and introduction, results and discussion are not fully separated.\nEven though in general, the paper is well-written, throughout the paper there are a some typos and grammar errors, e.g. (but not limited to):\nLine 87: ….only the three genes DFFB, …\nLine 98: locus = loci\nLine 441-442: Consequently, the origin ….. to the plasma. Incorrect sentence.\nLines 448-449: ….in which the cfDNA motifs…. Change into: …with the cfDNA motifs….\nLines 593-595: Recently, a newly….fetal genomes. Incorrect sentence.\nPlease go through the paper thoroughly to correct these."}, "author_response": "Responses: We sincerely thank the reviewer for the careful examination of our summary statistics and for pointing out this important issue. The INFO score shown in our previous figure represents the average imputation quality within each MAF interval, rather than the INFO score for individual variants. Although this was described in the Methods, it was not clearly emphasized in our previous response, which led to this misunderstanding.\nIn the revised version, we have explicitly added the INFO score as a variant-level filtering criterion. Specifically, in addition to the thresholds of MAF > 0.05 and HWE p > 1e-6, we now apply an INFO score > 0.4 cutoff in all GWAS analyses. Variants below any of these thresholds have been completely excluded. The updated GWAS summary statistics have been uploaded and are available at:\nhttps://drive.google.com/drive/folders/1Um6bPoBgw4Xokl677T0CDjD_8wxUjTvF?usp=sharing.\nFurthermore, we have re-performed all downstream analyses, including genetic correlation, Mendelian randomization, replication, and pathway-based analyses, based on the new GWAS results. The updated results have been incorporated into the revised manuscript. Importantly, the main findings and overall conclusions remain highly consistent with the original version. For example, the most significant associations continue to involve PANX1, DNASE1L3, and DFFB; replication rates have slightly improved; and the key enriched pathways remain those related to DNA fragmentation, caspase activity, and apoptosis. In addition, causal relationships inferred by MR analysis still remained between hematological phenotypes, particularly leukocyte traits, and cfDNA motif frequencies.\nTo facilitate comparison, we have provided a separate file “05-What was changed in this revision.docx” presenting the updated key results and figures alongside the previous versions, enabling a clearer visualization of the changes made.\nI address the reviewers' responses to my original comments below:\n1. Overall, the additional detail regarding imputation has adequately addressed my comments. In the response about INFO score thresholds, the authors provide a plot (S1C of their previous paper). It is unclear what the plotted INFO_SCORE is showing here; I'd assume median/mean(?) INFO score, per MAF bin. So while the average INFO score across all SNPs was >0.5, it isn't clear that each SNP had INFO > 0.5. Indeed, the summary statistics uploaded by the authors seemed to confirm that many variants had INFO < 0.5.\nResponses: We sincerely thank the reviewer for this comment. We apologize for the lack of clarity in our previous description. In the revised version, we have explicitly added the INFO score as an additional variant-level filtering metric. Specifically, we now require INFO > 0.4 for all retained variants, in addition to MAF > 0.05 and HWE p > 1e-6. All subsequent analyses, including GWAS, replication, Mendelian randomization, and genetic correlation, have been updated accordingly.\n2. Thank you for the additional clarification regarding trait transformation.\nResponses: Thank you.\n3. The additional explanation is appreciated.\nResponses: Thank you.\n4. The extra fields requested are available now in the uploaded summary stats but should also be included in the paper (e.g., table S1).\nResponses: We thank the reviewer for this comment. These extra fields of SNPs have now been included in Supplementary Table S3 of the revised manuscript.\n5. The additional information is welcome. Betas should be provided in table S3 in addition to p-values so consistency of direction of effect can be seen. It appears from Figures S8 and S9 that some variants have opposite directions of effect in the replication cohorts, compared to the discovery. If directions of effects are inconsistent, these shouldn't really be classified as a \"replicated\" association.\nResponses: We thank the reviewer for this comment. Accordingly, beyond assessing p-values, we have verified that the effect direction (beta) of the top SNP was consistent across discovery and replication studies for all replicated associations. The updated results have been incorporated and highlighted in the main text, and the beta values for these top SNPs in discovery and replication studies have been provided in Supplementary Table S3.\n6. The comparison of MR vs observed observations is interesting. Thank you for providing this; perhaps worth considering for a supplementary figure!\nResponses: We thank the reviewer for this suggestion. We have added this figure as Supplementary Figure S17.\n7. Perhaps worth trying to repeat relatedness estimation after filtering for poorly imputed variants.\nResponses: We thank the reviewer for this comment. We have recalculated the relatedness after filtering out variants with MAF, HWE, and INFO score. Below are the results:\nRelatedness Identical 1-degree 2-degree 3-degree Unrelated\n(pairs) twins relatives relatives relatives relatives\nGenotype 220 259 49,764,924 237,790,277 104,878,440\ndata without\nfiltering\nGenotype 220 207 8,748,364 161,135,145 222,550,184\ndata after\nfiltering\nThe new results showed that the numbers of predicted identical twins and 1-degree relatives barely changed, but the numbers of 2-degree and 3-degree relatives were dramatically reduced, and conversely, the number of unrelated relatives were significantly increased. The new results are more reasonable and improved, although an excess of related pairs still remains. We have added and highlighted the following sentence in the Discussion section:\nWhen applying additional quality control, using more reliable variants, such as those filtered by INFO score, results in more reasonable estimates of relatedness.\nReviewer #2 (Remarks to the Author):\nI am happy with the answers to my comments.\nJust a minor comment:\nIn you answer you state 'QUITL2, which allows simultaneous imputation of maternal and fetal genomes'.\nI believe this method is QUILT2.\nResponses: We thank the reviewer for the meticulousness. We have revised it as “QUILT2”\nReviewer #3 (Remarks to the Author):\nEven though the paper was rejected in the first place, the authors submitted a revised version. The authors did a tremendous piece of work, and thoroughly revised the paper. I would like to thank the authors for their efforts!\nAs far as my previous comments are concerned: the authors responded sufficiently to these. My first comment, however, is, despite the response of the authors, not completely resolved, but I leave it to the editor how to proceed with this. The paper is well-written, but it is still very lengthy. Moreover, the results section is full of background information and partly even discussion. I understand why the authors have chosen this format, as otherwise the rationale behind the different analyses would not be understood and the reader would be completely lost, but the fact is that this makes the paper very lengthy and introduction, results and discussion are not fully separated.\nResponses: We thank the reviewer for this valuable comment. We have carefully reviewed both the main text and the supplementary materials and have condensed the content accordingly. Specifically, we have removed repetitive descriptions across sections and transferred detailed methodological explanations to the Supplementary Information. Although the main text has been substantially shortened, the overall manuscript remains relatively long due to the comprehensive nature of the analyses and results.\nEven though in general, the paper is well-written, throughout the paper there are a some typos and grammar errors, e.g. (but not limited to):\nLine 87: ….only the three genes DFFB, …\nLine 98: locus = loci\nLine 441-442: Consequently, the origin ….. to the plasma. Incorrect sentence.\nLines 448-449: ….in which the cfDNA motifs…. Change into: …with the cfDNA motifs….\nLines 593-595: Recently, a newly….fetal genomes. Incorrect sentence.\nPlease go through the paper thoroughly to correct these.\nResponses: We thank the reviewer for this comment. We have carefully rechecked the entire manuscript to ensure that all typographical and formatting errors have been corrected."}, {"round": 3, "reviewer_comments": {"reviewer_1": "Thank you for addressing the concerns raised in my previous review. I'm pleased to see that the analyses have been re-done with the INFO filter, along with the updated MAF and HWE p-value filters. These changes, as anticipated, have indeed resulted in a more robust study with the increased replication rate. It's great to see that the majority of findings still hold. This revised version provides a much stronger proof of concept for the utility of cfDNA in genetic studies.\nI do have a couple of final minor points the authors should check:\nOn page 9, lines 240-41: The phrase \"among others\" following the list of replicated genes (PANX1, DFFB, and PADI4) seems superfluous and should be removed.\nOn page 10, lines 255-258 (and in the response to previous comment 5): The statement about consistent effect directions across studies isn't fully supported by Figure S8 and Table S3. There appear to be a small number of significant associations showing opposite effect directions in one of the replication cohorts (I think all are the natural population cohort) compared to the discovery study. This nuance should be accurately reflected in the text."}, "author_response": "Responses: We sincerely thank the reviewer for all the time and efforts dedicating to our work. With the constructive and insightful comments, our work has been significantly improved. We are deeply grateful for the support and recognition.\nI do have a couple of final minor points the authors should check:\n1. On page 9, lines 240-41: The phrase \"among others\" following the list of replicated genes (PANX1, DFFB, and PADI4) seems superfluous and should be removed.\nResponses: We thank the reviewer for this comment. We have revised and condensed the sentence as follows:\nOverall, three loci (PANX1, DFFB, and PADI4) were consistently replicated across all three studies, indicating shared genetic associations between pregnant and non-pregnant individuals.\nWe also noticed that the preceding paragraph (Lines 236–244) could be further condensed. The original version is as follows:\nIn the Chinese NIPT study, 14 out of 15 (93.3%) loci were replicated at a replication significance threshold of 1e-5, including PANX1, DFFB, DNASE1L3, and MSR1 (Figure 2C, Table S2). Among the 15 significant loci, the Dutch study successfully replicated 14 (93.3%), including PANX1, DFFB, DNASE1L3, and DNASE1L1 (Figure 2C, Table S2). Furthermore, 3 loci (20.0%), including PANX1, DFFB, and PADI4, were successfully replicated in the natural population study (Figure 2C, Table S2), demonstrating that the associations between cfDNA motif features and these genes are not limited to the pregnant population but are also present in the general population.\nWe have revised and highlighted this paragraph as follows:\nAt a replication significance threshold of 1e-5, 14 (93.3%), 14 (93.3%), and 3 (20.0%) loci were replicated in the Chinese NIPT, Dutch, and natural population cohorts, respectively (Figure 2C, Supplementary Data 2). Specifically, PANX1, DFFB, DNASE1L3, and MSR1 were replicated in the Chinese NIPT cohort; PANX1, DFFB, DNASE1L3, and DNASE1L1 in the Dutch cohort; and PANX1, DFFB, and PADI4 in the natural population cohort. The replication rates were notably higher in the pregnant cohorts than in the non-pregnant cohort, likely due to larger sample sizes. These findings suggest that the identified genetic associations are not limited to pregnant population but also in the general population.\n2. On page 10, lines 255-258 (and in the response to previous comment 5): The statement about consistent effect directions across studies isn't fully supported by Figure S8 and Table S3. There appear to be a small number of significant associations showing opposite effect directions in one of the replication cohorts (I think all are the natural population cohort) compared to the discovery study. This nuance should be accurately reflected in the text.\nResponses: We thank the reviewer for this valuable comment. We completely understand the confusion arising from our statement that “the effect directions of their top SNPs were consistently identical between the discovery study and the replication study”, as Figure S8 and Table S3 show a few points in the second and fourth quadrants. To clarify, this statement referred specifically to the top SNPs within each replicated locus. The points appearing in the second and fourth quadrants correspond either to non-replicated SNPs (hollow dots in Figure S8) or to non–top SNPs within replicated loci. For the 256 and 44 top SNPs (representing motif–locus associations) replicated in the Chinese NIPT and natural population cohorts, respectively, the effect directions were indeed consistent between the discovery and replication studies.\nTo avoid further confusion, we have revised Figure S8 to display only the top SNPs for the significant association in discovery study (solid dots indicate replicated, while hollow dots indicate non-replicated) and have bolded the beta values of the replicated motif-locus associations (top SNPs) in Table S3. We have also clarified the corresponding description in the main text as follows:\nFor the 256 and 44 replicated motif-locus associations in the Chinese NIPT and natural population cohorts, respectively, we further examined whether their top SNPs showed consistent effect directions between the discovery and replication studies (Supplementary Figure 8, Supplementary Data 3). The scatter plots demonstrated that all replicated motif–locus associations (solid dots) exhibited consistent effect directions across the corresponding studies."}]

Decision letter