Phenome-wide association study in 25,639 pregnant Chinese women reveals loci associated with maternal comorbidities and child health.

Conceptualization, Q.L. and Yulin Zhou; methodology, J.G., T.Z., Ying Zhou, and Q.L.; software, J.G., C.X., H.F., and Q.C.; formal analysis, J.G.; resources, Q.G., T.Z., Z.X., J.X., D.J., Y.Y., W.X., H.Z., and A.H.; data curation, J.G. and Q.G.; writing – original draft, J.G. and Q.L.; writing – review & editing, J.G., Q.G., T.Z., C.X., Z.X., H.F., Q.C., Ying Zhou, J.X., D.J., Y.Y., X.W., H.Z., X.J., Yulin Zhou, and Q.L.; visualization, J.G.; funding acquisition, J.G., Q.L., and Yulin Zhou.

The current study enrolled 48,734 pregnant women who gave birth between 2015 and 2020 at the Xiamen University-affiliated Maternity and Child Health Care Hospital ( Figures 1 and S1 ; Table 1 ). We genotyped 5,440,758 unique common germline variants (5,421,997 SNPs and 18,761 insertions or deletions [indels]) from NIPT ( STAR Methods ; Figure S1 ). The average calling rate of each variant was 9.43% (0.0144%–43.3%). We inferred the ancestry of the cohort by comparing it to the reference populations from the 1000 Genomes Project (1KGP). We observed that all 48,734 individuals share the same genetic background as CHX ( Figure S2 ). Furthermore, the genotypes’ allele frequencies strongly correlate with ChinaMAP 30 (ρ = 0.971, Figure S3 ) and CHX in the 1KGP (ρ = 0.969). In addition, we observed significantly higher genotyping concordance in sequences from mothers who undergo multiple NIPT (83.9%–88.1%) than those that have a single test performed (74.1%–78.2%, p  < 2.2 × 10 −16 ), suggesting that our genotyping method is reliable. To cope with the missing genotypes, we imputed a total of 37,482,377 germline variant genotypes and selected 5,957,600 variants (5,407,342 SNPs and 550,258 indels) based on the reference panel of the 1KGP variants with R 2  > 0.5 and minimum allele frequency > 0.01. Figure 1 A schematic view of this study In this study, we report an association study crossing genome and phenome in a population of pregnant women and neonates of Han Chinese based on ulcWGS-NIPT data and EHRs of the Health Commission of Xiamen. In addition, we use SAIGE for maternal PheWAS and perform further analysis based on the maternal pleiotropic TALs using a Bayesian model to identify the genetic variants simultaneously associated with maternal and newborn phenotypes. Table 1 Dataset characteristics NIPT dataset ( N  = 48,734): median (SD; IQR) Number of participants with ICD10 codes pregnant women ( n  = 25,639) children ( n  = 14,149) Gestation weeks of NIPT 16.57 (2.77; 14.6–18.9) – Age 33 (4.78; 29–36) 2.68 (0.468; 2–3) Gender female ( n  = 25,639) female ( n  = 6,593; 46.6%), male ( n  = 7,556; 53.4%) BMI (kg/m 2 ) 21.8 (3.15; 20.1–23.9) – Total fully specified ICD10 codes ( N  > 30) 268 133 SD, standard deviation; IQR, interquartile range; BMI, body mass index. A schematic view of this study In this study, we report an association study crossing genome and phenome in a population of pregnant women and neonates of Han Chinese based on ulcWGS-NIPT data and EHRs of the Health Commission of Xiamen. In addition, we use SAIGE for maternal PheWAS and perform further analysis based on the maternal pleiotropic TALs using a Bayesian model to identify the genetic variants simultaneously associated with maternal and newborn phenotypes. Dataset characteristics SD, standard deviation; IQR, interquartile range; BMI, body mass index. We also retrieved 38 germline whole-arm copy-number variations (CNVs) in the cohort, corresponding to a frequency of 3.90 × 10 −5 –3.74 × 10 −3 . As the sequences for NIPT were a mixture of maternal and fetal DNA, the fetal fraction ranged from 0.00584 to 0.401, also in the normal range 31 , 32 , 33 ( Figure S4 ). We obtained the disease phenotypes from the EHRs of the mothers and children, respectively ( STAR Methods ). After quality control, we included 268 maternal and 133 neonatal and early-childhood phenotypes with an overlap of 84 phenotypes, ranging in prevalence from 0.121%–33.4% in the mothers to 0.219%–66.4% in the children ( Tables 1 and S1 ). Of the maternal phenotypes, 32.1% were obstetric- (O00-O99) and genitourinary-related diseases (N00-N99). Of the children’s phenotypes, 27.8% were respiratory-related (J00-J99) and digestive-system-related diseases (K00-K93). We evaluated the associations between 5,957,600 genotyped common germline variants and 268 maternal disease phenotypes with pregnancy age and five top principal components of ancestry as covariates in 25,639 individuals ( STAR Methods ). To control for false discoveries, we downloaded 315,029 genetic associations (loci to traits) from the GWAS Catalog 34 and selected a set of 5,022 associations involving trait-associated SNPs encompassing the 4,826 genotyped SNPs and 76 phenotypes retrieved from our NIPT cohorts. We thus defined the 5,022 associations as “benchmarking associations,” which represent the known GWAS results in the variants included in the present study, including 151 associations exclusively reported in Chinese population. Our analysis revealed that the replication rates and significance levels (SLs) of the benchmarks increased with statistical power ( Figure 2 A). Notably, with a minimal statistical power of 0.8 and an SL of 0.05 ( STAR Methods ), our data reproduced 90.9% (20/22) of all benchmarking associations and 100% (2/2) of associations reported in Chinese population. Some of the associations included gestational diabetes mellitus (O24; 11:92965261:C:G, p  = 6.15 × 10 −10 , odds ratio [OR] = 1.12), obesity (E66; 16:53779455:T:G, p  = 5.59 × 10 −8 , OR = 1.26), asthma (J45, 10:9022753:T:C, p  = 1.59 × 10 −4 , OR = 1.37), thyrotoxicosis (E05, 5:77239986:T:G, p  = 2.00 × 10 −4 , OR = 0.826), psoriasis (L40, 6:31276305:C:T, p  = 3.59 × 10 −4 , OR = 4.95), esophagitis (K20, 4:76499631:C:T, p  = 5.65 × 10 −4 , OR = 1.46), endometriosis (N80, 7:49030398:G:T, p  = 6.37 × 10 −4 , OR = 1.95), and hypothyroidism (E03, 4:148710641:G:T, 8.52 × 10 −4 , OR = 0.716). Figure 2 Maternal PheWAS replication of GWAS Catalog SNP-phenotype associations (A) Each dot represents the -log10( p value) of a single SNP-phenotype association in maternal PheWAS. Different dot colors represent the different classes of maternal phenotypes. The colored squares and lines represent the replication rates of known GWAS associations ( p  < 0.05) in different populations at different statistical powers. (B) The boxplot shows the significance levels of the maternal TALs representing known GWAS association (benchmark set) increase with statistical power. The proxies used for the TALs were based on different thresholds of LD-R 2 (>0.3, >0.5, and >0.8). (C) The fold-of-enrichment of maternal PheWAS TALs in known GWAS associations of different populations at statistical power > 0.8 and LD-R 2  > 0.8. (D) The UpSet plot shows overlaps among maternal PheWAS TALs ( p   0.8. Maternal PheWAS replication of GWAS Catalog SNP-phenotype associations (A) Each dot represents the -log10( p value) of a single SNP-phenotype association in maternal PheWAS. Different dot colors represent the different classes of maternal phenotypes. The colored squares and lines represent the replication rates of known GWAS associations ( p  < 0.05) in different populations at different statistical powers. (B) The boxplot shows the significance levels of the maternal TALs representing known GWAS association (benchmark set) increase with statistical power. The proxies used for the TALs were based on different thresholds of LD-R 2 (>0.3, >0.5, and >0.8). (C) The fold-of-enrichment of maternal PheWAS TALs in known GWAS associations of different populations at statistical power > 0.8 and LD-R 2  > 0.8. (D) The UpSet plot shows overlaps among maternal PheWAS TALs ( p   0.8. The most significant tag SNPs in GWAS are typically located in the same linkage disequilibrium blocks (LD blocks) as the causal variants. Therefore, we expanded the benchmark SNPs to LD blocks defined by thresholding the squared coefficient of correlation (LD-R 2 ) in CHX of 1KGP. As a result, we found that the SLs of the LD blocks containing benchmark SNPs were highest at LD-R 2  > 0.8 ( Figure 2 B). Moreover, the SLs of the LD blocks also increased with the statistical power ( STAR Methods ). Therefore, we chose thresholds of statistical power (>0.8) and LD-R 2 (>0.8) to best replicate benchmark associations. Notably, with these criteria, our data exhibit enrichment of known GWAS associations, with the enrichments increasing with the SLs, from 0.885 (SL = 1 × 10 −3 ) to 2.01 (SL = 5 × 10 −8 ) in CHX ( Figure 2 C). Together, data from NIPT reproduced most known genetic associations and effectively controlled for confounding factors. After carefully benchmarking our findings with known genetic associations, we incorporated an additional filtering step to avoid genotyping errors in ulcWGS. Specifically, we selected the trait-associated-SNPs by absolute difference of pairwise R 2 (<0.01) between our cohort and the reference cohort (CHX of 1KGP; STAR Methods ). This approach ensured that the pairwise R 2 values of the resulting SNPs are consistent with the reference cohorts ( Figure S5 ). Finally, we yielded 2,883 SNPs associated with 26 maternal phenotypes at p  < 5 × 10 −8 and a statistical power greater than 0.8. Among these SNPs, 98.9% ( n  = 2,853) were non-coding, including 37.5% ( n  = 1,081) intergenic variants and 53.6% ( n  = 1,544) intronic variants. Only 1.1% ( n  = 30) of the SNPs were located in exons, affecting 21 genes ( Table S2 ). Furthermore, 99.5% ( n  = 2,869) were near known GWAS Catalog SNPs within a distance of less than 50 kb ( Figure S6 ). Subsequently, we aggregated these SNPs into 442 TALs based on LD-R 2 (>0.8). By comparing our resulting TALs with previously published association studies, we found that 351 TALs have not previously been reported in the GWAS Catalog, 34 GeneATLAS, 35 PheWAS, 36 and LabWAS 37 ( Figure 2 D; Table S2 ). These 351 TALs were associated with 23 disease categories ( Table S2 ). We also identified numerous pleiotropic TALs, which were associated with at least two phenotypes in our dataset (344 out of 442; Figure 3 A; Table S2 ). For instance, 10q11.22 ( n  = 13), 7q22.1 ( n  = 13), and 2q13 ( n  = 14) exhibited associations with at least ten maternal phenotypes. Additionally, the majority of these pleiotropic TALs were associated with diseases related to reproductive and childbirth conditions ( Figure 3 B). Of note, 96.6% of the pleiotropic TALs were associated with genitourinary conditions, such as female infertility (N97) linked to 1:147155826:C:G ( p  = 2.93 × 10 −32 , OR = 1.33) and 2:108999920:T:C ( p  = 7.65 × 10 −43 , OR = 1.39), female genital prolapse (N81) linked to 6:61238993:A:G ( p  = 4.55 × 10 −19 , OR = 0.634), and irregular menstruation (N92) linked to 2:108857756:G:C ( p  = 2.18 × 10 −20 , OR = 1.20); 58.6% of the pleiotropic TALs were linked to diseases occurring during pregnancy, childbirth, and the puerperium, such as miscarriage or stillbirth (O02) linked to 10:48576703:C:G ( p  = 4.03 × 10 −22 , OR = 1.42), spontaneous abortion (O03) linked to 2:108989311:G:A ( p  = 5.04 × 10 −20 , OR = 1.37), and “maternal care for other conditions predominantly related to pregnancy” (O26) linked to 10:48577013:G:T ( p  = 3.99 × 10 −21 , OR = 1.16); and 30.8% of the pleiotropic TALs were associated with endocrine, nutritional, and metabolic diseases, such as deficiency of other nutrient elements (E61) linked to 7:67431233:CA:C ( p  = 1.05 × 10 −12 , OR = 1.64) and polycystic ovary syndrome (PCOS; E28) linked to 2:109267108:T:G ( p  = 3.02 × 10 −11 , OR = 1.42). Figure 3 The associations of maternal PheWAS between TALs and different phenotypes (A) The Manhattan plot shows the genome-wide association between TALs and maternal phenotypes at statistical power > 0.8 and LD-R 2  > 0.8. Each dot represents a phenotype association at each TAL, while different dot colors represent different classifications of phenotypes. Finally, different dot sizes represent the odds ratio (OR). Numbers at the bottom indicate chromosomes, while the dashed line indicates the statistical significance ( p  = 5 × 10 −8 ). (B) The UpSet plot shows the number of unique TALs in a different classification of phenotypes at p   0.8, and LD-R 2  > 0.8. (C) The association count between each maternal phenotype and TAL at p   0.8, and LD-R 2  > 0.8. Notably, different colors represent the various categorizations annotated with published association studies. “Replicated” refers to variants at LD-R 2  > 0.8 associated with the same trait annotated in known databases, “related” refers to variants affecting the genes known to associate with the trait of interest in Harmonizome, “reported” refers to variants in LD block (LD-R 2  > 0.8) with known TALs but to different phenotypes, and “novel” refers to variants with an association not previously reported. (D–H) The LocusZoom plot displays association p values on the -log10 scale on the vertical axis. At the same time, the chromosomal position is indicated along the horizontal axis for (D) female infertility (N97), (E) thalassemia (D56), (F) salpingitis and oophoritis (N70), (G) polycystic ovary syndrome (E28), and (H) maternal care for intrauterine death (O36). The associations of maternal PheWAS between TALs and different phenotypes (A) The Manhattan plot shows the genome-wide association between TALs and maternal phenotypes at statistical power > 0.8 and LD-R 2  > 0.8. Each dot represents a phenotype association at each TAL, while different dot colors represent different classifications of phenotypes. Finally, different dot sizes represent the odds ratio (OR). Numbers at the bottom indicate chromosomes, while the dashed line indicates the statistical significance ( p  = 5 × 10 −8 ). (B) The UpSet plot shows the number of unique TALs in a different classification of phenotypes at p   0.8, and LD-R 2  > 0.8. (C) The association count between each maternal phenotype and TAL at p   0.8, and LD-R 2  > 0.8. Notably, different colors represent the various categorizations annotated with published association studies. “Replicated” refers to variants at LD-R 2  > 0.8 associated with the same trait annotated in known databases, “related” refers to variants affecting the genes known to associate with the trait of interest in Harmonizome, “reported” refers to variants in LD block (LD-R 2  > 0.8) with known TALs but to different phenotypes, and “novel” refers to variants with an association not previously reported. (D–H) The LocusZoom plot displays association p values on the -log10 scale on the vertical axis. At the same time, the chromosomal position is indicated along the horizontal axis for (D) female infertility (N97), (E) thalassemia (D56), (F) salpingitis and oophoritis (N70), (G) polycystic ovary syndrome (E28), and (H) maternal care for intrauterine death (O36). We compared our TALs with previously published association studies. Among the significant associations reported, we identified two associations labeled as "replicated" ( Figure 3 C; Table S2 ; STAR Methods ). For instance, 1:107809190:T:C (1p13.3) was associated with hypothyroidism (E02; p  = 2.98 × 10 −8 , OR = 1.19) and 11:92959391:C:T (11q14.3) was associated with diabetes mellitus in pregnancy (O24; p  = 2.87 × 10 −9 , OR = 1.12). Furthermore, we classified an additional 242 associations (involving 225 TALs) as "related." Specifically, in this class of associations, the TALs were located in the known pathogenic genes of the associated traits. For example, 2:111056516:T:C (2q13) is associated with female infertility (N97; p  = 1.29 × 10 −35 , OR = 1.35) and located in ACOXL ( Figure 3 D), and 16:268762:G:A (16p13.3), associated with thalassemia (D56; p  = 6.12 × 10 −23 , OR = 3.21), is located in FAM234A ( Figure 3 E). Notably, most of the affected phenotypes in this class ( n  = 233; 96.3%) were related to female infertility (N97) and irregular menstruation (N92). Moreover, 35.0% of the related TALs (105/225) were known expression quantitative trait loci (eQTLs) of corresponding genes, based on the eQTL Catalogue 38 Finally, we identified 371 associations (85 TALs) as “reported,” where the TALs were known to be associated with different phenotypes. For example, 2:108897145:A:G (2q13) was associated with salpingitis and oophoritis (N70, p  = 4.16 × 10 −12 , OR = 1.62; Figure 3 F), which is an exonic variant of RANBP2 , a pathogenic gene of ovarian diseases and cancer. These findings suggest that many known TALs are potentially pleiotropic. Of note, in addition to the other three classes of associations, we identified 1,273 novel associations (344 TALs associated with 23 disease categories), which were discovered in our dataset ( Table S2 ). The majority of these novel associations were related to diseases in the genitourinary system and those occurring during pregnancy. Moreover, some of these novel associations were independently validated by independent studies. For instance, the intronic variant 10:48551720:G:A (10q11.22) in ARHGAP22 has been associated with PCOS (E28; p  = 9.42 × 10 −9 , OR = 1.40; Figure 3 G), and ARHGAP22 has been recognized as a crucial gene in PCOS. 39 Additionally, TRRAP is implicated in early embryonic lethality and defects in cell cycle progression. 40 We discovered that the intronic variant, 7:98909792:C:CCTTT (7q22.1), of TRRAP was associated with maternal care for intrauterine death (O36, p  = 3.30 × 10 −8 , OR = 1.21; Figure 3 H). Furthermore, NPHP1 has been implicated in male infertility. 41 , 42 Then, we found four intronic variants of NPHP1 associated with female infertility (N97), such as 2:110160461:A:AT (2q13, p  = 4.87 × 10 −20 , OR = 1.24; Figure 3 D). These findings shed light on the genetic underpinnings of various maternal phenotypes and underscore the importance of these novel associations in understanding complex diseases and their potential implications for female health. In addition, we observed that many of the maternal diseases showed a polygenic nature ( Figure 3 C; Table S2 ). For instance, female infertility (N97) is associated with 402 unique TALs. In comparison, female genital prolapse (N81) is associated with 353 unique TALs. These findings support the notion that the cumulative effects of multiple genetic variants are often linked to complications during pregnancy, childbirth, and puerperium. 43 In addition to SNPs and indels, we evaluated germline whole-arm CNV events, comprising 13 microduplications and eight microdeletions, to investigate their association with 268 maternal phenomes ( Figure S7 ). Our analysis revealed 13 microduplications and three microdeletions that were significantly associated with maternal diseases ( p  < 0.05; Figure 4 ; Table S3 ). Notably, microduplications in 18p ( p  = 7.31 × 10 −16 , OR = 61.74), 18q ( p  = 9.40 × 10 −24 , OR = 97.4), and 21q ( p  = 2.24 × 10 −33 , OR = 88.1) were found to be associated with “maternal care for known or suspected fetal abnormality and damage” (O35). Additionally, a microduplication in 13q was linked to disorders of the lacrimal system (H04; p  = 2.94 × 10 −5 , OR = 90.2). Of note, the associations for 13q have been previously reported. 44 Regarding the microdeletions, we found that a microdeletion in 19q was associated with retained placenta (O73; p  = 2.13 × 10 −5 , OR = 128.6), one in 14q was associated with biliary tract diseases (K83; p  = 3.32 × 10 −5 , OR = 91.5), and another in 22q was associated with thalassemia (D56; p  = 3.99 × 10 −5 , OR = 36.0). Figure 4 The associations between germline CNVs and maternal phenotypes The association between germline CNVs and maternal phenotypes. Each dot represents a CNV-phenotype association, and different dot colors represent a different classification of maternal phenotypes. Finally, different dot sizes represent the OR. Numbers to the left indicate the significant level: upward for microduplications and downward for microdeletions. Labels at the top indicate whole-chromosome arms. The associations between germline CNVs and maternal phenotypes The association between germline CNVs and maternal phenotypes. Each dot represents a CNV-phenotype association, and different dot colors represent a different classification of maternal phenotypes. Finally, different dot sizes represent the OR. Numbers to the left indicate the significant level: upward for microduplications and downward for microdeletions. Labels at the top indicate whole-chromosome arms. Maternal-neonatal comorbidities are known from many studies. 45 , 46 , 47 Thus, we investigated whether there is an association between maternal pleiotropic TALs and complications in children. Using a Bayesian model, we evaluated the associations between 91,826 maternal TALs ( p   0.8) and 58 neonatal phenotypes in 14,149 individuals. The model estimated the ratios of allele frequencies (α) between the case and control groups ( STAR Methods ). We identified 21 maternal pleiotropy SNPs in 16 TALs associated with 35 neonatal phenotypes ( Figure 5 A; p   0.8; Table S4 ). Notably, we reported two TALs associated with the same complications in mothers and children. 19:28958117:T:C (19q12) was associated with dermatitis (L30) in mothers ( p  = 9.93 × 10 −6 , OR = 1.24) and children ( p  = 0.049, the ratio of reference allele frequencies between case and control: α = 7.12). Similarly, 15:57500986:A:G (15q21.3) was associated with acute sinusitis (J01) in mothers ( p  = 3.12 × 10 −4 , OR = 1.57) and children ( p  = 0.046, α = 7.31). Furthermore, we noticed that some maternal pleiotropic TALs were also pleiotropic in children ( Figure 5 A). For instance, 3:45186585:G:A is associated with maternal hypothyroidism (E03, p  = 3.21 × 10 −5 , OR = 0.703), while in children, it is associated with childhood gastroenteritis and colitis (K52, p  = 0.0265, α = 9.20), gastroenteritis and colitis of unknown origin (A09, p  = 0.0377, α = 7.99) and acute tonsillitis (J03, p  = 0.0373, α = 8.02). Similarly, 8:94799830:T:C was associated with maternal diseases of complicating pregnancy, childbirth, and puerperium (O99, p  = 1.91 × 10 −5 , OR = 0.752), which is also associated with childhood pityriasis alba (L30, p  = 0.0295, α = 8.82), melaena (K92, p  = 0.0398, α = 7.82), and acute laryngotracheitis (J04, p  = 0.0406, α = 7.75) in children. Figure 5 The associations between maternal pleiotropic TALs ( p  < 1 × 10 −3 ) and neonatal phenotypes (A) The heatmap shows the association between maternal and neonatal phenotypes that shared the same TAL. The color indicates significant associations between maternal pleiotropic TALs and neonatal phenotypes. The asterisk indicates that children’s disease prevalence was significantly different (Fisher test p  < 0.05) in different maternal genotypes and disease statuses. (B–F) Bar plot shows the disease prevalence of each neonatal phenotype at different maternal genotypes and disease statuses. Labels at the top indicate the maternal phenotypes and the Fisher test p value, while labels at the bottom indicate the genotypes of the specific SNP. The numbers on the left indicate the prevalence of a specific phenotype in children. Different colors represent the different disease statuses, and the numbers labeled on the histogram represent the number of individuals. The associations between maternal pleiotropic TALs ( p  < 1 × 10 −3 ) and neonatal phenotypes (A) The heatmap shows the association between maternal and neonatal phenotypes that shared the same TAL. The color indicates significant associations between maternal pleiotropic TALs and neonatal phenotypes. The asterisk indicates that children’s disease prevalence was significantly different (Fisher test p  < 0.05) in different maternal genotypes and disease statuses. (B–F) Bar plot shows the disease prevalence of each neonatal phenotype at different maternal genotypes and disease statuses. Labels at the top indicate the maternal phenotypes and the Fisher test p value, while labels at the bottom indicate the genotypes of the specific SNP. The numbers on the left indicate the prevalence of a specific phenotype in children. Different colors represent the different disease statuses, and the numbers labeled on the histogram represent the number of individuals. We further verified the association between maternal pleiotropic TALs and neonatal traits based on the observed disease prevalence in each genotype ( Figures 5 and S8 ). Our data revealed that common pediatric respiratory diseases were associated with eight unique maternal TALs. For example, hypothyroidism (E03) showed an association with pulmonary diseases in adults, 48 particularly a higher risk of pneumonia (J18). 49 , 50 We found that maternal TALs (3p21.31) of hypothyroidism (E03) were also linked to neonatal pneumonia (J18; α = 7.39; Figure 5 B), as well as acute neonatal bronchitis (J20; α = 7.39; Figure 5 C). Moreover, atopic dermatitis (L30) was associated with increased risk of systemic infections 51 and bronchitis in early childhood. 52 We observed that acute neonatal bronchitis (J20) was also associated with maternal TALs of atopic dermatitis (L30, 3p25.2, α = 8.36; Figure 5 D). Furthermore, our study showed that neonatal pulmonary diseases were associated with maternal diseases complicating pregnancy. For example, neonatal pneumonia (J18) was associated with “maternal care for other conditions predominantly related to pregnancy” (O26, 7q22.1, α = 7.44; Figure 5 E), and neonatal bronchitis (J40) was associated with maternal TALs of “maternal diseases classifiable elsewhere but complicating pregnancy, childbirth, and puerperium” (O99, 8q22.1, α = 7.35; Figure 5 F). These findings indicate potential associations between maternal and neonatal diseases, offering valuable insights into shared genetic factors that may contribute to comorbidities in both populations, especially in pediatric respiratory diseases.

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Qiyuan Li ( [email protected] ). This study did not generate new unique reagents. • The raw genotype data have been deposited at Genome Variation Map and are publicly available as of the date of publication. Accession numbers are listed in the key resources table . Phenotype data reported in this paper will be shared by the lead contact upon request. • All original codes reported in this paper are publicly available on GitHub: https://github.com/xmbd/PheWAS-across-Chinese-pregnant-women-from-ulcWGS and also has been deposited at Zenodo: https://zenodo.org/records/11365274 . DOIs are listed in the key resources table . • Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. The raw genotype data have been deposited at Genome Variation Map and are publicly available as of the date of publication. Accession numbers are listed in the key resources table . Phenotype data reported in this paper will be shared by the lead contact upon request. All original codes reported in this paper are publicly available on GitHub: https://github.com/xmbd/PheWAS-across-Chinese-pregnant-women-from-ulcWGS and also has been deposited at Zenodo: https://zenodo.org/records/11365274 . DOIs are listed in the key resources table . Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Indeed, even with the recent establishment of consortia of large genetic and phenotypic databases, genome- or phenome-level association studies are still costly. NIPT based on ulcWGS is the most widely utilized prenatal genetic test. 22 Its widespread adoption globally provides an excellent opportunity for conducting association studies on large populations of pregnant women and newborns. In this study, we presented a comprehensive association study crossing whole-genome and phenome data, leveraging genetic information from NIPT and a regional EHR database. Our research provides insights into the genetic and phenotypic association landscape within a CHX population of pregnant women and their newborns. It is worth noting that despite significantly compromised coverage and fetal effects, we successfully replicated 86.4% (317/367) of known associations in the GWAS Catalog at a power level above 0.3. These results demonstrate the potential of utilizing NIPT data to validate and enhance associations identified through GWAS and PheWAS. Furthermore, our data shed light on the genetic factors underlying maternal and neonatal diseases. This study contributes to a deeper understanding of the genetic basis of various health conditions in the context of pregnancy and childbirth. GWAS can characterize the genetic variants associated with a single trait at a whole-genome level. At the same time, PheWAS aim to reveal the associations between a set of genetic variants and many traits. Our analysis provides a comprehensive view of the association landscape between genetic variants and multiple traits related to non-gestational and gestational complications, as well as maternal-neonatal morbidities. Interestingly, our data reveal that a significant proportion, 82.9% (922/1,112), of the unique TALs (LD-R 2  > 0.8) exhibit pleiotropic effects, influencing multiple disease traits. It is worth noting that many known TALs were found to be associated with disease traits not tested initially in previous studies. For example, the TAL associated with red blood cells (16:572565:A:G) was associated with thalassemia (D56, p  = 2.09 × 10 −9 , OR = 0.734), while the TAL of fasting blood glucose (11:92935660:C:T) was associated with gestational diabetes mellitus (O24 p  = 3.70 × 10 −8 , OR = 1.11). These findings are classified as reported in our study and can offer valuable insights into the biological background of the TALs. Nevertheless, the related class associations provide additional evidence to the known related genes ( Figures 3 D–3H). Notably, many genetic variants associated with the same diseases also serve as eQTLs for the related genes ( Table S2 ). Such findings further reinforce the relevance and significance of these genetic associations in understanding disease mechanisms. Maternal and newborn health are intricately connected, 53 , 54 and our study has provided important insights into the genetic basis of many maternal-neonatal comorbidities. Firstly, the discovery of these genetic associations can inform future better prevention strategies for pregnant women and children. Our study reported potential risk loci for maternal-neonatal comorbidities. For instance, we observed an association between a non-coding variant, 19:28958117:T:C (19q12), and atopic diseases in both mothers ( p  = 1.06 × 10 −5 , OR = 1.24) and children ( p  = 0.049, α = 7.12), which is one of the major causes of disease burden in this population worldwide. 55 , 56 Pneumonia (J18) is known to have a high incidence among children (22.1% in our study), and we identified an association between maternal pleiotropic TALs and pneumonia. For example, the intronic variant 7:102215622:G:A in the CUX1 was associated with maternal care for abnormality of pelvic organs (O34; p  = 1.10 × 10 −8 , OR = 1.15), hemorrhage in early pregnancy (O20; p  = 4.79 × 10 −5 , OR = 1.07), and female genital prolapse (N81; p  = 8.97 × 10 −10 , OR = 0.735). Our findings enable healthcare professionals to better diagnose and mitigate risks in both mothers and their children. Furthermore, the risk loci reported also help to better understanding the genetic risk in this population and improve clinical genetic counseling programs in obstetrics and pediatrics. Then, the present study demonstrated the potential of circulating DNA sequences as a powerful tool for diagnosis and genetic screening. Currently, NIPT is only used for detection of trisomy, while based on the results of this study, we can envision more extensive clinical usage of maternal circulating DNA in the prenatal diagnoses of genetic disorder and rare diseases. Finally, our study also provided a scalable analytical pipeline for association studies based on EHRs and low-coverage circulating DNA sequences, which can be applied to larger populational research. In addition, the biological background of the reported TALs remains unexplored. Investigation into the underlying biological processes will offer great opportunity to identify novel therapeutic targets and unlock new avenues for personalized medicine in maternal and neonatal healthcare. While the current study has made valuable discoveries, it is essential to consider the following limitations when interpreting the results. Firstly, the read length and coverage of NIPT sequences resulted in minimal calling rates for known variants, leading to a constrained sample size for testing single variants. While essential for controlling false positive associations, this strict quality control filtering approach may have also excluded some true novel associations that did not meet the filtering criteria. Nevertheless, the analysis still managed to uncover numerous new associations, thanks to the utilization of EHR-based phenome data and the sheer size of the population. Secondly, it is essential to note that the study exclusively focused on pregnant women, which inevitably introduced a bias toward obstetrics and genitourinary complications. Consequently, the data may only partially reveal genetic associations corresponding to uncommon disease traits that are underrepresented in the cohort. Thirdly, we specifically focused on disease diagnoses during pregnancy, resulting in highly unbalanced case and control samples for many phenotypes. To address this challenge, we applied the SAIGE method to correct for the inflation resulting from the unbalanced case-control ratio. Lastly, it should be highlighted that the NIPT sequences also contain information on somatic mutations in the population. Exploring and investigating the implications of these somatic mutations could be a promising direction for future studies in the field. In summary, we conducted a comprehensive association analysis at both the genome and phenome levels, leveraging NIPT and EHR data. We demonstrate the importance of accumulating NIPT data to inform GWAS and PheWAS while uncovering crucial potential genetic determinants related to maternal diseases and maternal-newborn comorbidities. These findings have the potential to complement existing GWAS and PheWAS efforts. They may facilitate the development of more targeted and personalized approaches to prenatal care and pediatric preventive strategies, ultimately improving maternal and neonatal health outcomes.

Genome- and phenome-wide association studies (GWAS and PheWAS) are both common and important methods used to identify the genetic determinants of human diseases. 1 GWAS, empowered by high-throughput sequencing technology, can detect linkages between many variants and single traits, while PheWAS identify pleiotropic loci simultaneously associated with multiple traits. The recent establishment of extra-large biospecimens collections, such as the UK Biobank and Global Biobank Meta-analysis Initiative, 2 , 3 can enable the genetic association studies of various diseases with statistical power. 4 In particular, mining electronic health records (EHRs) with “big data” analyses greatly expands the spectrum of phenotypes 5 , 6 , 7 and enhances the capacity to uncover associations between diseases and genetic variations, advancing our understanding of the genetic background of diseases. 8 , 9 Nevertheless, challenges in EHR-based PheWAS include limited case numbers and imbalanced data, which may inflate type I errors and impact the sensitivity and specificity of the results. 10 , 11 Additionally, with the expansion of the test space, PheWAS requires intensive computing power. 12 To address these challenges, new methods, such as SAIGE and GMMAT, 10 , 13 , 14 have been proposed to analyze the imbalanced samples due to different prevalence of diseases. Ultra-low-coverage whole-genome sequencing (ulcWGS; <1×) has been used in population genetic studies encompassing high genome coverage and low cost. 15 , 16 Non-invasive prenatal testing (NIPT), which uses ulcWGS, is a leading technology in prenatal screening for fetal chromosomal abnormalities. 17 NIPT has higher sensitivity and specificity than traditional maternal serum screening in large-scale clinical trials. 18 , 19 Recently, it has been also applied to the diagnosis of genetic disorders such as sex chromosome aneuploidies. 20 , 21 NIPT has become one of the most successful clinical applications of circulating fetal cell-free DNA-based tests. 22 In addition, the massive amount of DNA sequences generated by NIPT can provide population-wise genomic information for association studies. 23 , 24 , 25 The genetic information obtained via NIPT can also be utilized to identify loci associated with maternal or fetal phenotypes, such as cancer. 26 , 27 , 28 However, a major obstacle is NIPT’s ultra-low sequencing depth (around 0.0561–0.371×) and short read length (35 or 28 bp), leading to limited mapping quality. Although large consortia of genetic data make high-quality GWAS and PheWAS feasible, the high cost of the data acquisition process must be considered before conducting full-scale association studies. 29 So far, few association studies have focused on maternal diseases and maternal-newborn comorbidities, especially in the Chinese population. 3 Meanwhile, NIPT data with matched, well-structured EHRs can help reveal genetic variants associated with specific disease phenotypes and complement findings on maternal diseases and maternal-newborn comorbidities obtained through GWAS and PheWAS 25 analyses. Here, we report an association study of both the genome and phenome in a Han Chinese (CHX) population of pregnant women and newborns, with NIPT data and EHRs collected by the Health Commission of Xiamen. The study covered 5,957,600 genetic variants and 317 recorded international codes for diseases (ICD10). We report 2,883 trait-associated-SNPs, including 442 trait-associated loci (TALs) associated with 26 pleiotropic loci. A further analysis of the maternal pleiotropic SNPs using a Bayesian model identified 21 genetic variants simultaneously associated with maternal and newborn phenotypes. Our findings provide a comprehensive view of the genetic background of non-gestational and gestational complications. We also identify germline variants underlying the maternal-neonatal association of morbidities. These findings help fill the blank in the phenome-wide genetic determinants in CHX women, inform future GWAS and PheWAS, and improve maternal and infant health.

The authors declare no competing interests.

REAGENT or RESOURCE SOURCE IDENTIFIER Deposited data BioSample This paper BioProject: PRJCA017926 Raw variants and genotypes data This paper GVM: GVM000558 Human reference genome NCBI build 37, GRCh37 Genome Reference Consortium http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/human/ Human reference genome NCBI build 38, GRCh38 Genome Reference Consortium http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/human/ 1000 Genomes Project Phase 3 data 1000G Consortium 2015 Nature ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502 GWAS catalog Sollis, E. et al. 34 https://www.ebi.ac.uk/gwas/downloads GeneATLAS Canela-Xandri, O. et al. 35 http://geneatlas.roslin.ed.ac.uk/ PheWAS catalog Denny, J.C et al. 36 https://phewascatalog.org/ LabWAS catalog Goldstein, J.A. et al. 37 https://phewascatalog.org/labwas Harmonizome Rouillard, A.D. et al. 57 https://maayanlab.cloud/Harmonizome/ eQTL Catalog Nurlan Kerimov et al. 38 https://www.ebi.ac.uk/eqtl/ ChinaMAP Cao Y. et al. 30 http://www.mbiobank.com/ Genome Aggregation Database (gnomAD) Karczewski, K.J. et al. 58 https://gnomad.broadinstitute.org/ ICD10 codes and their related phecodes This paper; Zenodo Data Table S1 .tsv; Zenodo: https://zenodo.org/records/11365274 Maternal trait-associated-SNPs This paper; Zenodo Data Table S2 .tsv; Zenodo: https://zenodo.org/records/11365274 Maternal trait-associated-CNVs This paper; Zenodo Data Table S3 .tsv; Zenodo: https://zenodo.org/records/11365274 Maternal trait-associated-SNPs associated with neonatal diseases This paper; Zenodo Data Table S4 .tsv; Zenodo: https://zenodo.org/records/11365274 Maternal clinical information This paper; Zenodo Data mother_clinical_info.tsv; Zenodo: https://zenodo.org/records/11365274 Neonatal clinical information This paper; Zenodo Data child_clinical_info.tsv; Zenodo: https://zenodo.org/records/11365274 Software and algorithms BWA version 0.7.17-r1188 H. Li et al. 59 https://sourceforge.net/projects/bio-bwa/files/ GATK version 4.1.4.1 Auwera, G.A. et al. 60 https://github.com/broadinstitute/gatk/releases Picard version 1.119 Broad Institute https://broadinstitute.github.io/picard/ Bvcftools version 1.9 Danecek P et al. 61 https://samtools.github.io/bcftools/howtos/install.html CrossMap version 0.5.4 Zhao, H. et al. 62 https://crossmap.sourceforge.net/ Plink2 Chang, C.C et al. 63 https://www.cog-genomics.org/plink/2.0/ Eagle2 Loh et al. 64 http://data.broadinstitute.org/alkesgroup/Eagle/downloads/ Minimac4 Das, S. et al. 65 https://genome.sph.umich.edu/wiki/Minimac4 WisecondorX Raman, L. et al. 66 https://github.com/CenterForMedicalGeneticsGhent/WisecondorX/ PREFACE Raman, L. et al. 67 https://github.com/CenterForMedicalGeneticsGhent/PREFACE SAIGE Zhou, W. et al. 14 https://saigegit.github.io/SAIGE-doc/ R version 3.6.2 R Software https://www.r-project.org/ Python version 3.7.4 Python Software Foundation https://www.python.org R packages genpwr Moore, C.M. et al. 68 https://cran.r-project.org/web/packages/genpwr/index.html Bayesian model for Maternal Pleiotropic Loci associate with neonatal phenotypes This paper 08_PheWAS_fetal.md; Zenodo: https://zenodo.org/records/11365274 The study included 48,734 women who underwent NIPT between 12 and 31 weeks of pregnancy. Women were registered at XMUMCH between 2015 and 2020. All the pregnant women or their authorized relatives were fully informed of the study and provided written consent. This study was approved by the ethics committee/institutional review board (IRB) of XMUMCH. The enrolled subjects’ phenome data (diagnoses in ICD10) were retrieved from the EHR system managed by the HCXM, which included outpatient and inpatient information from various hospitals in Xiamen. Owing to privacy and data access restrictions, we were confined to using only the primary ICD10 diagnosis codes for analysis. To elucidate the relationship between ICD10 codes and phecodes, we have detailed the correspondence between these codes and their associated phecodes, providing comprehensive information for each in Table S1 . We used the disease category code to define the phenotypes and filtered out the categories of non-heritable diseases, such as external causes of morbidity and mortality; symptoms, signs, and abnormal clinical and laboratory findings not elsewhere classified; and injury, poisoning, and additional consequences of external causes. Due to privacy and data access restrictions, we were limited to use the primary ICD10 diagnosis codes. For the mothers, we retrieved both the outpatient and inpatient EHR ranging from two years pre- and post-partum. In addition, we retrieved all the EHR (outpatient and inpatient) for the children up to three years after birth. Cases and controls were named based on the occurrence of a given disease category to define binary traits for the association study. Finally, disease categories with either low case or control numbers ( N  < 30) were also excluded for insufficient statistical power. The raw and incremental data were cleaned and standardized to form a Medical Data Asset. The corresponding datasets were saved in the Privacy-Preserving format. Details of the sequencing protocol were described in a previous study. 24 The cleaned reads of 31,063 samples were trimmed to 35bp. Next, we aligned the reads against the human genome ref. 38 (hg38) with bwa using the single-end read alignment option. 59 Potential PCR duplicates were removed with Picard (v1.119). Finally, we recalibrated the base quality with GATK4. 60 The other 17,671 NIPT data were trimmed to 28bp read length and aligned against hg19. We converted coordinates from hg19 to hg38 using CrossMap. 62 We used BCFtools to call and normalize the variants controlling for base quality ≥20 and mapping quality ≥ 30. 61 We filtered the variants with a mappability uniqueness score of less than one. 24 We then determined the genotype based on the allele’s depth as follows: (1) reference allele depth (RD) > 0 and alternate allele depth (AD) = 0, homozygote of the reference allele coded as 0; (2) RD > 0 and AD > 0, heterozygote coded as 1; (3) RD = 0 and AD > 0, homozygote of the alternate allele coded as 2. We annotated the germline variants using ANNOVAR and selected the common germline variants with minimum allele frequency (MAF) > 0.01 in Chinese population from ChinaMAP, 30 and East Asian population (EAS) from The Genome Aggregation Database (gnomAD). 58 We retrieved germline variants from 48,734 NIPT data. We compared the MAF of the cohort with those of the Chinese population in ChinaMAP, 30 as well as the East Asian population (EAS) in The Genome Aggregation Database (gnomAD). 58 We imputed genotypes based on the 1KGP reference panel using Eagle2 64 and Minimac4 65 with inclusion criteria of MAF>0.01 and imputation score R 2  > 0.5. 69 To verify their accuracy, we explored genotype concordance rates between NIPT the same mother who undertake multiple NIPT and that from single tests. As for CNVs, we first selected a group of healthy subjects with pregnancy ages under 35 and no recorded history of present illness upon NIPT. Of note, these patients had no ICD10 records and no history of present illness (e.g., median- or high-risk screening for Down syndrome, spontaneous abortions, and abnormal B-ultrasound). Next, we used WisecondorX 66 to identify copy-number alterations for 31,063 individuals with 646 reference samples, while the remaining 17,671 individuals with 683 reference samples had bin sizes of 100kb. Finally, we identified the whole-arm CNVs of which fragment lengths were at least 80% of the chromosome arms. We inferred the ancestry of 48,734 NIPT individuals using principal component analysis (PCA) based on the genotyping of common autosomal variants and reference populations, specifically Han Chinese (CHX, n  = 208), Dai Chinese (CDX, n  = 109), and Japanese (JPT, n  = 105) from 1KGP. To improve accuracy and save computing power, we selected the variants located on autosomal chromosomes with MAF >0.1, calling number among the top 1%, and the observed allele frequency (OAF) within the 95% confidence interval of the alternate allele frequency using the binomial distribution. We performed population classification using Support Vector Machine (SVM) based on the top two principal components (2-PCs). We inferred the fetal DNA fraction of chromosome Y (FFY) of the male fetus and defined the FFY of the female fetus as 0 70 using the following equation. We performed PREFACE to build a fetal fraction (FF) statistic model based on CNVs and FFY. 67 (Equation 1) F F Y = % c h r Y − f e m a l e % c h r Y m a l e % c h r Y − f e m a l e % c h r Y Where % c h r Y refers to the fraction of bases mapped to Y chromosome, f e m a l e % c h r Y and m a l e % c h r Y refer to the median % c h r Y in female and male fetuses, respectively. In the current analysis, m a l e % c h r Y was set to 0.00173 according to previous studies. 71 We performed PheWAS for all maternal variants against all phenotypes with pregnancy age and top-5 PCs of population structure as a covariate. We used SAIGE to adjust for unbalanced case-control ratios and sample relatedness in the data. 14 To validate our PheWAS, we downloaded 315,029 genetic associations (loci to traits) from GWAS Catalog. 34 And we defined the associations involve trait-associated-SNPs which were retrieved from our NIPT cohorts “benchmarking associations”, which represent the known GWAS results in the variants included in the present study. We also calculated the power for each pair of variants and phenotypes using the R package "genpwr", 68 based on sample size (N), case rate, observed minimum allele frequency (OMAF), odds ratio (OR), and the underlying genetic model. The significance level was set to 0.05. 36 Test results with low statistical power (0.01 and imputation score R 2  > 0.5, n  = 5,957,600) using plink2 63 in our cohort and the Chinese Han population ( n  = 208) of 1KGP. Then, to further control for false discoveries caused by genotyping errors, we compared the distributions of pairwise-R 2 between our cohort and the reference cohort (CHX of 1KGP), and selected the SNP-pairs of which the absolute difference of pairwise-R 2 between the two cohorts are less than 0.01, thus yielded 2,358,854 SNPs. Then, we used Pearson’s correlation coefficient to evaluate the consistency of LD-structures of the filtered SNP set in our cohort and the reference cohort, and we kept the intersection of previous significant SNPs associated with maternal traits and the filtered SNP set ( n  = 2,358,854) for further analyses. All the significant variants from the PheWAS results were manually reviewed and classified into four categories, as follows: (1) “Replicated” refers to variants in LD-R 2 (>0.8) with variants associated with the same trait annotated in public GWAS databases; (2) “Related” refers to the tag-SNP (i.e., the SNP with highest significance) or the proxies in LD-R 2 (>0.8) are located in the genomic region of the corresponding pathogenic gene of the same disease as documented in Harmonizome 57 ; (3) “Reported” refers to variants in LD-block (LD-R 2 >0.8) with known TAL but to different phenotypes; (4) “Novel” refers to significant associations between trait-associated-SNPs and phenotypes that were uniquely discovered in our dataset and have not been reported elsewhere. Furthermore, these SNPs their proxies in LD (LD-R 2  > 0.8) have not been reported to be associated with any phenotypes in existing databases, including GWAS Catalog, 34 GeneATLAS, 35 PheWAS 36 and LabWAS. 37 We also performed PheWAS for all maternal germline whole-arm CNV against all phenotypes with pregnancy age and top-5 PCs of population structure. We derived a Bayesian model to evaluate the association between variants and children’s phenotypes from mixture sequence data such as NIPT with a very low fraction of fetal DNA. RD and AD refer to the reference and alternate allele depth in the raw genotyping data before imputation. We divided the children into cases (positive diagnosis) and controls (negative diagnosis) for each phenotype. We assumed that the RDs and ADs followed beta-binomial distribution in each cohort. (Equation 2) ( R D , A D ) ( 1 − δ 0 0 1 + δ ) ( 2 ( 1 − ϕ ) 0 0 2 ϕ ) ∼ B B ( R D , R D + A D , P , θ ) (Equation 3) P = p + π ∗ μ (Equation 4) p ∈ A F C H X ± 3 ∗ σ (Equation 5) σ = A F C H X ∗ ( 1 − A F C H X ) c a l l i n g n u m b e r (Equation 6) α = P c a s e P c o n t r o l (Equation 7) log ( α ) ∼ N ( 0 , τ 2 ) Where: δ refers to base-calling error (0< δ <1 × 10 −3 , ϕ refers to alignment error rate (0.4< ϕ <0.6); P refers to the observed allele frequency measured from mixture sequences (0< P <1); p refers to the expectation of reference maternal allele frequency, of which the boundaries are defined by to allele-frequency in CHX ( A F C H X ) as eq. d.; σ refers to the standard deviation; π refers to the fetal DNA fraction; μ refers to the mutation rate in the children; α refers to the ratio of reference allele frequencies between case and control, following a log-normal distribution of N ( 0 , τ 2 ) .

The authors would like to thank Professor Jiahuai Han, Dr. Huanming Yang, and Dr. Bangchang Xie for critical reading of and feedback on our work. We thank the Health Commission of Xiamen and BaseBit ( https://www.BaseBit.Ai ) for support in phenome data application and processing. We thank BGI ( https://www.bgi.com ) for providing sequencing data. J.G. was supported by a grant from the 10.13039/501100001809 National Natural Science Foundation of China (award number 82203420 ). Q.L. was supported by the 10.13039/501100012166 National Key Research and Development Program of China (award numbers 2022YFA1303200 and 2022YFA1303201 ). Yulin Zhou was supported by Major 10.13039/100017357 Science and Technology Projects of Fujian Province (award number 2021YZ034011 ) and major scientific research projects of 10.13039/501100014125 Fujian Provincial Health Commission (award number 2021ZD01002 ).