This two-stage case-control genome-wide association study (GWAS) involved unrelated Chinese participants recruited from both West China Hospital (hospital-based population cohort) and affiliated community clinics (natural population cohort) through advertisements, word of mouth, and referrals from clinicians or other researchers. Advertisements offered the opportunity for PFO screening at no cost and did not indicate monetary compensation. There was no specified other purpose or population selection for study participants in order to maximize subject recruitment. Participants were initially referred to the Department of Cardiology at West China Hospital for clinical interviews and PFO screening. Baseline characteristics, including age, gender, ethnicity, educational level, body mass index (BMI), smoking, alcohol consumption, regular coffee intake, and medical history, were recorded during the clinical interview. Inclusion criteria were: 1) aged 4–80 years; 2) understand and cooperate with contrast transthoracic echocardiography (cTTE). Exclusion criteria were: 1) self-reported history of severe organic diseases, or tumor; 2) non-Han ethnicity to minimize population stratification; 3) presence of any heart disease other than PFO.

There were 581 participants with PFO and 2,646 participants without PFO recruited in discovery cohort between 10 March 2020 and 18 January 2023, and 2,924 samples passed filtering (see later). Due to imbalances in baseline characteristics between cases and controls in the discovery cohort (Table 1), propensity score matching (PSM) was performed to avoid bias caused by imbalanced sample number and baseline characteristics. Propensity scores were calculated using covariates including age, gender, educational level, BMI, smoking, alcohol consumption, and coffee intake. Cases and controls were matched 1:1 using the nearest-neighbor method via the MatchIt package. Finally, 517 cases and 517 controls were used in the first discovery stage. There were 263 individuals recruited in discovery cohort between 7 February 2023 and 19 March 2024, and no samples with more than 10% missing genotypes or gender nonconformity.

Written informed consent was obtained from all participants, and the study was approved by the Ethics Committee of West China Hospital of Sichuan University (No. 2020145) and registered at the Chinese Clinical Trial Register (ChiCTR1900024623) in accordance with Helsinki Declarations.

Transthoracic echocardiography (TTE) was initially performed using 1–5 MHz or 3–8 MHz multiplane transducers and a Philips IE33 ultrasound device to detect any other cardiac diseases. After TTE showed no structural or functional abnormalities, cTTE was performed. Specifically, a contrast medium containing microbubbles (an agitated solution generated by 1 mL of blood, 1 mL of air, and 8 mL of saline) was injected into the antecubital veins (Silvestry et al., 2015). Participants were assessed at rest, during a Valsalva maneuver (expiratory pressure at 60 mm Hg measured by manometer), and while coughing. As described in previous studies, enough Valsalva maneuver could enhance PFO detection (Zhao et al., 2019). Following Chinese expert consensus on the diagnosis of PFO(Chinese Consensus Expert Group On The Clinical Application Of Transesophageal Echocardiography, 2022), we measured the expiratory pressure during the Valsalva test by a standard manometer. Patients were asked to exhale into a mouthpiece connected with the manometer trying to achieve 60 mmHg and maintain the expiratory pressure for 15 s. Then patients were instructed to release on command after arrival of contrast in the right atrium. Subjects were allowed to practice the Valsalva maneuver prior to the formal test. PFO was diagnosed if three or more detected microbubbles per frame appeared in the left heart within three cardiac cycles (Maffe et al., 2010).

For the size of the PFO was quantified based on the number of detected microbubbles per frame in the left heart: Grade 1, less than 10 microbubbles; Grade 2, 10–30 microbubbles; Grade 3, more than 30 microbubbles (Maffe et al., 2010).

Peripheral blood samples (5 mL) were collected from the basilic vein into standard EDTA-containing collection tube. DNA was extracted from whole blood using a MagPure Blood DNA LQ kit with a KingFisher Flex automated nucleic acid purification system (Thermo Fisher Scientific, Waltham, MA, United States) following the manufacturer’s instructions. A total of 500–1,000 ng of DNA per sample was used for DNA preparation. The integrity of gDNA was determined through gel electrophoresis. DNA were fragmented on Covaris LE220-Plus. PCR products were purified using VAHTS™ DNA Clean Beads (Vazyme Biotech, Nanjing, China) in a Tecan Evo instrument. Quality control was performed using the Agilent Bioanalyzer 4200 (Agilent Technologies, Santa Clara, CA, United States) to ensure an average fragment size of 460 bp. DNA concentrations were determined with a Tecan Infinite 200 PRO and were required to be more than 5 ng/μL. DNB preps of clinical samples generated in Tecan 480 were sequenced on the ultra-high-throughput DNBSEQ-T7 platform (MGI, Shenzhen, China).

Whole-genome sequencing reads were mapped to the human reference genome (GRCh38/hg38). The GATK best practice pipeline was applied to call variants and indels (McKenna et al., 2010). The PLINK program (v 1.9), KING (v 2.3.0) and R statistics (v 4.1.0) were used for quality control procedures. Samples were removed if they exhibited: (1) gender inconsistencies (n = 5 in cases, n = 22 in controls), (2) genotype call rates less than 90% or outlying heterozygosity (n = 41 in cases, n = 55 in controls), or (3) third-degree relative or closer relationships (n = 18 in cases, n = 162 in controls). SNPs were filtered if: (1) the call rate was less than 90%, (2) the minor allele frequency (MAF) was less than 1%, or (3) significant deviation from Hardy-Weinberg equilibrium (p < 10⁻⁴). Variants were annotated using the ensemble variant effect predictor (VEP) (McLaren et al., 2016). All annotation processes were conducted based on GRCh38 genome coordinates.

Amplification reactions were performed in a total volume of 15 μL, containing 2 μL 10 × buffer (Roche), 1 μL dNTP, 0.5 μL forward primer (10 μM), 0.5 μL reverse primer (10 μM), 0.15 μL Taq DNA polymerase (Roche), 1 μL genomic DNA and 9.85 μL pure water. PCR conditions were: 2 min at 95°C, 30 cycles of 30 s at 95°C, 30 s at 65°C, and 30 s at 72°C, with a final extension for 2 min at 72°C. Alternatively, PCR was performed using PrimeSTAR Max DNA Polymerase (Takara) with polymerase activation at 98°C for 30 s, followed by 35 cycles at 98°C for 10 s and 67°C for 60 s. Primers for mutation confirmation with Sanger sequencing are shown in the Supplementary Materials. PCR fragments were purified from agarose gels, and sequencing was performed using the BigDye Terminator cycle sequencing kit (Applied Biosystems). DNA fragment analysis was performed on an ABI 3730 Genetic Analyzer (Applied Biosystems) according to the manufacturer’s instructions. Sequencing profiles shown in this manuscript were generated by SnapGene Viewer 4.0 (Chicago, IL, United States).

Two public datasets from articles by Asp et al. (2019) and Cui et al. (2019) containing accessible single-cell data of human fetal hearts were included. For each dataset, gene expression UMI count values were log-normalized using Seurat v4 NormalizeData with a scale factor of 10,000. After centering and scaling with Seurat ScaleData, principal component analysis (PCA) was performed on variable genes with Seurat RunPCA and embedded in two-dimensional Uniform Manifold Approximation and Projection (UMAP) plots with Seurat RunUMAP on 10 (data from Asp et al. (2019)) or 30 (data from Cui et al. (2019)) principal components. Clusters showing expression levels of CNOT2, KCNMB4, MLLT10, FRG1 and IGBP1 were visually selected. The SingleR package was utilized to identify the predominant cell types. Considering PFO closure is a process of endothelial-to-mesenchymal transition mediated fibrosis (Elliott et al., 2014), pseudotime analysis was performed on mesenchymal stem cells using to predict possible differentiation trajectory using the Monocle 2 package. Genes of interest were modeled as smooth, nonlinear function of pseudotime, and gene expression changes along this pseudotime were tested.

Student’s t test and χ² test were used to compare continuous outcomes and categorical outcomes between cases and controls in both discovery and validation cohorts. All tests were two-tailed and p < 0.05 was considered significant. Continuous data are presented as mean ± standard deviation (SD) and categorical data as numbers and percentages (%). All statistical analyses were conducted using R version 4.1.0 unless otherwise specified.

PLINK was used to perform genotypic association analyses. For binary traits (without PFO vs. with PFO), logistic regression was used with adjustment for the top three principal components. A p value threshold of less than 1 × 10⁻⁶ in the discovery cohort was used for suggestive genome-wide association (Schaid et al., 2018) to select candidate SNPs for validation follow-up. Quantile–quantile (QQ) plots were generated to evaluate the extent to which the observed GWAS p values deviated from the null hypothesis. To estimate heritability attributable to genome-wide assayed genetic variation, the genome-wide complex trait analysis (GCTA) package was used (Yang et al., 2011). A linear regression model for the size of PFO was further assessed with PLINK in the discovery cohort. The GTEx eQTL Calculator was used for expression quantitative trait loci analysis.

The study flow from initial participant screening through recruitment and selection is summarized in Figure 1. After matching demographic covariates, the discovery cohort finally included 517 PFO cases and 517 demographically matched controls (Table 1). PCA showed no stratification among participants (Supplementary Figure 1). The Manhattan plot for the main GWAS analysis is presented in Figure 2A, and the dash line marks the genome-wide suggestive significance threshold (p < 1 × 10⁻⁶). There were 3 SNPs (rs78526049, rs1313755017, rs1227675732) reaching genome-wide significance (p < 5 × 10⁻⁸) (Table 2). Among 11 SNPs reaching suggestive significance, eight were intergenic variants (rs1313755017 and rs1227675732, rs79366272 and rs78526049, rs62206791 and rs62206790 were close with a distance not more than 200bp), and three were intron variants (Table 2). The QQ plot (Figure 2B) shows the distribution of observed and expected p values, with minimal genomic control inflation (λ = 0.998). Basic heritability estimates are presented in Supplementary Table 1. After GWAS, all identified variants were confirmed through manual inspection of bam files in Integrative Genomics Viewer (IGV). Following this verification and the removal of false positives (rs62142007), Sanger sequencing was performed (Supplementary Figure 2). It was confirmed that rs1227675732, rs62206790, rs879176184, rs13115019 and rs57922961 were true positives. According to the results of GWAS in discovery cohort, DUX4L35 (rs1227675732) (OR = 2.903, 95% CI 1.961–4.297, p = 3.05 × 10⁻⁸) showed the most significant among all true positive SNPs. A second independent association was observed at rs62206790 (OR = 2.780, 95% CI 1.864–4.146, p = 2.02 × 10⁻⁷) in MLLT10P1. SNPs in FRG1HP (rs879176184) (OR = 2.724, 95% CI 1.822–4.073, p = 4.30 × 10⁻⁷) and IGBP1P5 (rs13115019) (OR = 2.437, 95% CI 1.702–3.488, p = 5.80 × 10⁻⁷) were also increased risk for PFO. Conversely, the rs57922961 (T > C) variant showed association with a reduced risk of PFO (OR = 0.508, 95% CI 0.388–0.666, p = 6.82 × 10⁻⁷) (Table 2).

Among the 517 subjects with PFO, 121, 125 and 271 cases were categorized as Grade 1, 2 and 3 respectively. We observed that mutations in DUX4L35 (rs1227675732), MLLT10P1 (rs62206790), FRG1HP (rs879176184), IGBP1P5 (rs13115019) mainly presented in PFOs with large shunt (Grade 3) (Supplementary Figure 3A). The SNP rs57922961, associated with a reduced risk of PFO, was mainly found in PFOs with moderate-to-large shunts (Grade 2–3) (Supplementary Figure 3A). Exploratory analyses of quantitative measures of PFO were also conducted. None of the SNPs analyzed met conventional criteria for genome-wide significance (Supplementary Figure 3B). However, rs62206790, rs1227675732, rs13115019 and rs879176184 showed nominal evidence of association in PFO size (p < 0.01) (Supplementary Table 2). The QQ plot illustrated observed p values were smaller than those expected, likely due to insufficient sample size for detecting effects on PFO size (Supplementary Figure 3C).

Replication genotyping was performed in a validation cohort of 111 PFO cases and 152 controls (Table 3). The flow chart of patient enrollment of validation cohort is shown in Figure 3A. There were 66.7% cases (74/111) with large shunt. The variations of rs62206790, rs879176184, rs13115019 were more prevalent in participants with PFO, and the variation of rs57922961 was less common in participants with PFO (Figure 3B). The trends of these four SNPs variation in PFO were recapitulated in the validation cohort, although the differences were not significant at the nominal level (p = 0.0576, p = 0.3972, p = 0.7383, p = 0.3883, respectively). However, individuals with ≥2 of these four site variations were much more likely to have PFO in the validation cohort (p = 0.0109). The rs1227675732 variation was not observed in the validation cohort.

Given that the four identified SNPs were in intergenic or intronic regions, their functional potential was annotated using publicly available databases. We focused on rs57922961 and its associated gene expression (Figure 4A), while MLLT10P1 (rs62206790), FRG1HP (rs879176184) and IGBP1P5 (rs13115019) could not be tested for eQTLs in whole blood or heart tissues according to GTEx Project datasets. This SNP was found to modify expression levels of KCNMB4, a gene 11.3 kb downstream of CNOT2, based on the GTEx database (Figure 4A). Figure 4B showed the expression levels of genes correlated with rs57922961 in different tissues. We found that the rs57922961-T allele was associated with higher expression levels in whole blood (p = 2.36 × 10⁻⁷) but was not significantly correlated with mRNA expression levels in heart atrium tissue (p = 0.565) (Figure 4C).

To further explore the function of the genes in heart development, we leveraged publicly available human single-cell RNA-seq data of heart tissue. The single-cell RNA sequencing data from Asp et al. (2019) were generated from 19 heart tissue sections from the 4.5 to 9 post-conception weeks, and the data from Cui et al. (2019) were generated from 18 human embryos’ hearts ranging from 5 to 25 weeks of gestation. Dimensionality reduction and cell type categorization were performed based on marker gene expression for the filtered cells in the two datasets separately (Figures 5A, B). As FRG1HP, MLLT10P1 and IGBP1P5 are pseudogenes, we replaced them with their homologous functional genes. Figures 5C, D revealed the expression of FRG1, MLLT10, IGBP1, CNOT2 and KCNMB4 in distinct populations of cells in 4.5–9-week and 5–25-week hearts, separately. Monocle2 pseudotime trajectory analyses were used to develop a developmental course of mesenchymal stem cells (MSC). Two branch points were identified in 4.5–9-week hearts (Figure 5E) with KCNMB4 and IGBP1 expression declining after an initial increase (Figure 5F). CNOT2, MLLT10 and FRG1 showed a drop in expression with a slight increase in the middle across pseudotime in 4.5–9-week hearts (Figure 5F). The same pseudotime analysis in the 5–25-week hearts also identified cell fates in an unsupervised manner (Figure 5G). It also identified heightened expression of IGBP1 at the beginning, with decreased expression along the pseudotime trajectory (Figure 5H). Restricting the analyses to hearts in the 5–25 W stage, similar results were observed with CNOT2, MLLT10, FRG1 and KCNMB4 showing slight increases during the descent process along pseudotime.