This study (Approval No.: BJH-2024004) was approved by the Ethics Committee of Beijing Jiaen Hospital. All participants provided informed consent, with procedures complying with the Declaration of Helsinki (1964) and amendments.

In February 2024, a couple with PREMBL was enrolled at Beijing Jiaen Hospital. Family history was documented and the couple completed clinical examinations, laboratory tests, and basal hormone analysis (Cobas® 6000, Roche Diagnostics).

Karyotyping (550-band G-banding), chromosomal microarray analysis (CMA, Cytoscan 750k), and whole-exome sequencing (WES, Illumina NovaSeq 6000) were performed on the couple.

Whole-exome sequencing (WES) was performed on the proband’s DNA sample to detect sequence variants, following a methodology previously described (Yang et al., 2019). Target regions were enriched using the Agilent Sure Select Human Exon Sequence Capture Kit (Agilent Technologies, United States). Enrichment efficiency of the DNA library was assessed by quantitative PCR. Library fragment size, distribution, and concentration were determined using the Agilent Bioanalyzer 2100 (Agilent Technologies, United States). Sequencing was performed on the NovaSeq6000 (Illumina, Inc.) platform with paired-end reads of approximately 150 bp. Libraries (normalized to ∼300 pM per sample) were sequenced using the NovaSeq Reagent kit. Raw sequencing reads (quality level Q30% > 90%) were aligned to the human reference genome (accession no: hg19/GRCh37) using the Burrows-Wheeler Aligner (BWA). PCR duplicates were marked and removed using Picard tools (v1.57). Variant calling was performed using the Verita Trekker® Variant Detection System (v2.0; Berry Genomics, China) in conjunction with the Genome Analysis Tool Kit (https://software.broadinstitute.org/gatk/). Variant annotation and interpretation were conducted following the joint guidelines of the American College of Medical Genetics and Genomics (ACMG) (Richards et al., 2015), utilizing ANNOVAR v2.0 (Wang et al., 2010) and the Enliven® Variant Annotation Interpretation System (Berry Genomics). Variant pathogenicity assessment was supplemented using data from the three frequency databases (1000G_2015aug_eas, https://www.internationalgenome.org; ExAC_EAS, http://exac.broadinstitute.org; gnomAD_exome_EAS, http://gnomad.broadinstitute.org), as well as the Human Gene Mutation Database (HGMD) pro v2019. Additionally, the Revel score (for pathogenicity prediction) (Ioannidis et al., 2016) and the pLI score (representing the tolerance for truncating variants) were applied.

Sanger sequencing, performed on an Applied Biosystems 3500DX Genetic Analyzer (Thermo Fisher Scientific, United States), was used to validate the identified variants.

The impact of missense variants on amino acid evolutionary conservation was analyzed using MEGA7 software (http://www.megasoftware.net/previousVersions.php; default parameters). The conservation score of specific amino acid residue positions was visualized using the WEBLOGO tool (http://weblogo.berkeley.edu/logo.cgi).

Homology modeling was performed using SWISS-MODEL, based on the crystal structure template PDB: AF-Q6TGC4-F1 (https://swissmodel.expasy.org), to generate three-dimensional (3D) structural models of the wild-type (WT) PADI6 protein and the mutant proteins (p.L237Afs*24 and p.E670Gfs*48).

Subsequently, molecular dynamics analysis (MD) were conducted on the wild-type (WT) and mutant PADI6 models (p.L237Afs*24 and p.E670Gfs*48) generated by Modeller 9v17 (Sali et al., 1995). Hydrogen atoms and N-/C-terminal patches were added using the CHARMM22 force field (MacKerell et al., 1998). The models were solvated in a TIP3P water box for hydration and charge neutralization, with the minimum distance between the model and the box boundary set to 13 Å. Simulations were performed using the NAMD 2.9 software under periodic boundary conditions (PBC) (Phillips et al., 2020). Key simulation parameters included: temperature maintained at 300 K, pressure at 1 atmosphere (atm), and time steps of 2 femtoseconds (fs). Electrostatic interactions were modeled using the particle mesh Ewald (PME) method, with a van der Waals interaction cutoff distance of 12 Å. The models underwent a three-stage pre-equilibration process (total duration: 600 ps), the last snapshots were selected as the beginning structures for 40 ns production simulations under unrestrained conditions.

To construct expression plasmids encoding wild-type (WT) or mutant PADI6 (p.L237Afs*24 and p.E670Gfs*48), total RNA was extracted from the proband’s mRNA sample. WT and mutant PADI6 cDNA sequences were amplified by RT-PCR. The resulting cDNA fragments were subcloned into the pET vector. After verification by Sanger sequencing, the fragments were cloned into the pcDNA3.1(+)-3×Flag vector. The final constructs were designated as PADI6-WT, PADI6-L237Afs*24 and PADI6-E670Gfs*48.

HEK 293T cells (ATCC, Cat# CRL-3216) were cultured and seeded in 24-well plates. Upon reaching appropriate confluency, cells were transfected with the respective plasmids using Lipofectamine 3000(Invitrogen, United States). Cells were harvested at 36 or 48 h post-transfection for subsequent analysis.

Cells were harvested 48 h post-transfection. Total RNA was extracted using the RNeasy Mini Kit (QIAGEN), followed by cDNA synthesis with the PrimeScript RT Reagent Kit (including gDNA Eraser, Takara). Gene expression levels were quantified using gene-specific primers, normalized to the reference gene GAPDH (MIM:138400). Real-time qPCR was performed on a 7500 Fast Real-Time PCR System (Applied Biosystems) with three technical replicates per sample. Data are presented as mean ± standard deviation (SD) from three independent experiments. The PCR cycle was as follows: 95 °C for 10 min, 1 cycle; 95 °C for 10 s + 60 °C for 30 s (fluorescence acquisition), 55 cycles. The expression values of each gene were normalized to the expression level of GAPDH using the 2^−ΔΔCT method.

At 36 h post-transfection, cells were fixed with 4% paraformaldehyde for 30 min, followed by permeabilization with 1% Triton X-100 at room temperature for 1 h. Samples were blocked in blocking buffer (3% BSA and 5% goat serum in PBS) for 1 h. Primary incubation was performed with anti-Flag mouse monoclonal antibody (1:400 dilution; Sigma-Aldrich, United States) at 4 °C overnight. After PBS washes, samples were incubated with Alexa Fluor-conjugated goat anti-mouse secondary antibody (1:400 dilution; Invitrogen, Carlsbad, CA) for 2 h at room temperature protected from light. Nuclei were counterstained with DAPI (Invitrogen, Carlsbad, CA) for 5 min. Coverslips were mounted using Fluoromount aqueous mounting medium (Sigma-Aldrich, St Louis, WA) and imaged under a laser scanning confocal microscope (Zeiss, Gottingen, Germany) with a ×63 oil immersion objective.

Cells were harvested 48 h post-transfection and lysed with 2× SDS lysis buffer containing 1 mM phenylmethanesulfonyl fluoride (PMSF) (Sigma, St Louis, WA), 0.2 mM β-mercaptoethanol, and protease inhibitor cocktail (Sigma, St Louis, WA). Equal protein lysates were resolved by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE), transferred onto polyvinylidene fluoride (PVDF) membranes, and subjected to immunoblotting with anti-Flag mouse monoclonal antibody. β-actin served as the loading control.

Statistical analyses were performed using Prism 8 software (GraphPad Software, San Diego, CA). Intergroup differences were evaluated by Student's t-tests. Statistical significance was defined as p ≤ 0.05.

The pedigree chart constructed through systematic collection of family history is presented in Figure 3A. The proband (II-1, female) and her spouse (II-2, male) were both 33 years old with a 5-year history of infertility. Body mass index (BMI) and basal hormone levels were within normal ranges for both individuals. Gynecological examination revealed right unilateral hydrosalpinx in II-1. For II-2, no erectile or ejaculatory dysfunction was reported, and semen kinetic analysis demonstrated normal sperm motility and morphological parameters. Both individuals denied a history of sexually transmitted infections (STIs), long-term medication use, or exposure to harmful chemicals or radiation.

The couple underwent two in vitro fertilization (IVF) cycles, both of which employed a GnRH antagonist protocol. In the first cycle, twelve oocytes were retrieved, eight of which were successfully fertilized and progressed to cleavage. The second cycle yielded fifteen oocytes, with nine achieving successful fertilization and undergoing cleavage. Crucially, in both cycles, all embryos that successfully fertilized and reached the cleavage stage exhibited developmental arrest on culture day 3 (at the 4 to 8-cell stage) (Figure 1).

Chromosomal analysis at 550-band resolution revealed a normal male karyotype (46,XY) (Figure 2A). The female karyotype was reported as 46,XX,?der(14)ins(14;?)(p12;?), suggesting an insertion of unknown origin at chromosome 14p12. Subsequent high-resolution chromosomal microarray analysis (CMA) confirmed a normal female genomic copy number profile arr(1-22,X)×2, leading to the classification of this karyotypic abnormality as a clinically insignificant benign polymorphism (Figures 2B,C).

WES identified no pathogenic variants in the male partner (Figure 3B). In the female proband, compound heterozygous frameshift variants in the PADI6 gene (RefSeq transcript NM_207421.4) were detected: c.707dupT (L237Afs*24) and c.2009_2010delAG (E670Gfs*48). Both variants induce a frameshift. The c.707dupT variant, inherited from the proband’s father, results in a premature termination codon (PTC) and predicted protein truncation (Figures 3A,I–1). Conversely, the c.2009_2010delAG variant, inherited from the proband’s mother, leads to a delayed termination codon and predicted protein elongation (Figures 3A,I–2). Following the American College of Medical Genetics and Genomics (ACMG) guidelines, these variants were classified as “Likely Pathogenic” (PVS1+PM2+PP4) and “Pathogenic” (PVS1+PM2+PP4), respectively. This classification was based on their loss-of-function mechanism, absence in population databases (gnomAD), and their phenotypic correlation with Preimplantation Embryonic Lethality 2 (PREMBL2; OMIM: #616814), an autosomal recessive disorder. Database searches confirmed that this specific compound heterozygous combination has not been previously reported.

Multiple sequence alignment of PADI6 protein sequences from seven species (Figure 4B). The results demonstrated that Leu237 (L237) and Glu670 (E670) exhibited 100% conservation scores across all analyzed species. Visualization of the ±12 amino acid regions flanking these residues using WEBLOGO revealed single dominant peak characteristics for both positions, indicating strong purifying selection during evolution.

Homology modeling via the SWISS-MODEL server was employed to generate three-dimensional structures of wild-type (WT) PADI6 and its variants p.L237Afs*24 (L237A) and p.E670Gfs*48 (E670G) (Figure 4C). The p.L237Afs*24 mutation localizes to the intermediate domain of PADI6, where the introduction of a premature termination codon results in a truncated protein backbone of 259 amino acids. This represents a 62.7% reduction in length compared to the 694-amino-acid WT protein and leads to complete ablation of the C-terminal domain (Figure 4A).

The p.E670Gfs*48 variant resides within the C-terminal domain (Figure 4A). This mutation extends the protein backbone to 716 amino acids, incorporating 22 additional residues at the C-terminus relative to WT. Structural analysis reveals the disruption of two critical hydrogen bonds maintained between Glu670 and other residues in the wild-type (WT) conformation. Consequently, this alteration induces localized structural differences at the C-terminal region.

Molecular dynamics analysis (MD) during 0–40 ns revealed that the L237A mutant exhibited unstable overall conformations (Figure 5A). The E670G variant maintained relatively stable global conformations but developed localized structural changes at the C-terminus during later simulation stages, which was corroborated by root mean square deviation (RMSD) trajectories and root mean square fluctuation (RMSF) analysis (Figures 5B–D).

Further RMSF analysis indicated enhanced flexibility of key residues and structural instability in critical local regions of L237A, specifically manifested by near-complete α-helix disappearance, reduced β-sheet content, and increased random coils/turns beyond residue 170 (Figure 5F). While E670G showed dynamic distributions similar to WT overall, it exhibited higher residue flexibility in certain regions. Pronounced structural fluctuations emerged beyond residue 660, characterized by destabilized α-helices and diminished β-sheets (Figure 5G).

Both variants demonstrated reduced hydrogen bonding between mutation sites and neighboring residues compared to WT (Figure 5E). In summary, L237A and E670G affect protein structural integrity and stability to varying degrees.

Statistical analysis of mRNA expression levels in the wild-type group (PADI6-WT), L237A mutant group (PADI6-L237Afs*24), and E670G mutant group (PADI6-E670Gfs*48) revealed a reduction in relative mRNA expression for both mutants, with concomitant significant impairment in the expression of mRNAs encoding PADI6-interacting proteins (Figure 6A).

Subsequent immunofluorescence experiments demonstrated that neither the L237A nor E670G variants altered the subcellular localization of the protein, with both mutants exhibiting cytoplasmic distribution identical to WT (Figure 6B). However, comparative analysis of relative fluorescence intensity revealed significantly attenuated fluorescence signals in both mutants (Figure 6D).

Further Western blot analysis detected bands for both WT and E670G mutants at approximately 80 kDa, whereas the L237A mutant band at around 30 kDa, consistent with pre-experimental molecular weight expectations (WT: 80 kDa; L237A: 31 kDa; E670G: 82 kDa) (Figure 6C). Concurrently, the band signal intensities across three groups corroborated the immunofluorescence results. Statistical quantification of relative protein expression levels further confirmed significantly diminished expression in both mutants (Figure 6E).

Collectively, based on the results of protein expression, molecular weight, abundance and interacting proteins, it was confirmed that both mutations affect PADI6 function.