Clinical manifestations

The general data, visual acuity results, fundus photographs, and FFA of eight affected individuals, including the proband and first-degree relatives from four families were shown in Table 1. The probands consisted of three males and one female, with a mean age of 12.25 years (range, 7 to 17 years), while the affected family members included one male and three females, with a mean age of 43.25 years (range, 42 to 45 years). BCVA ranged from finger count before the eye to 1.0. Notably, fundus manifestations varied significantly among individuals, even within the same family. 3 out of 4 probands were categorized as stage 3, while their parents were classified as stage 2. Specifically, the FFA of the probands from three families exhibited RD in either of the eyes, except for the proband from family 1, whose fundus only showed dilated blood vessels and fluorescein leakage. In contrast, their parents consistently showed no perfusion area, increased vascular branching, and stained vascular walls (Fig. 2).

Table 1 Summary of general information and phenotype characteristics of affected individuals from 4 FEVR familiesFig. 2

Fundus Photography and Fluorescein Angiography of Individuals Affected by FVER. Fundus photographs and fluorescein angiography images of probands in four different families are presented as increased vascular branches, fluorescein leakage, no perfusion area scleral buckling or paraoptic disc retinal atrophy, retinal detachment, while their parents only showed no perfusion area and increased vascular branches or fluorescein leakage. Abbreviations: FVER stands for familial exudative vitreoretinopathy, OD for the right eye, and OS for the left eye

Variants detected by targeted genetic sequencing

Four variants were detected in 4 out of 9 (44%) families (Table 2, Fig. 3A). Among these, two novel heterozygous variants (c.4289delC: p.Pro1431Argfs*8 and c.2073G > T: p.Trp691Cys) were identified in LRP5, while the other two variants (c.1801G > A: p.Gly601Arg in LRP5 and c.633 T > A: p.Tyr211* in TSPAN12) had been previously reported [50, 51].

Table 2 Computational analysis of variants in four probands with FEVRFig. 3

Genetic findings in four families and conservative analysis of a novel missense variant. A Pedigrees of the families with causative variants and the corresponding Sanger sequencing results of a heterozygous variant in the proband and his/her first-degree non-affected relative. Squares represent males, circles represent females, arrows indicate the proband, open symbols represent unaffected individuals, solid symbols represent affected individuals, ± represents heterozygous variant, and + / + represents wildtype. B Alignment of LRP5 homologous protein sequences from 26 species, using human LRP5 as a reference, reveals a remarkable conservation of the tryptophan residue at position 691 across these species. C The phylogenetic trees represent the evolutionary relationships of various species

The newly identified variant c.4289delC: p.Pro1431Argfs*8 in LRP5 caused a Pro to Arg substitution at position 1431, followed by premature termination at position 1437 due to eight frameshift amino acids. This resulted in a truncated protein product. The second novel variant, c.2073G > T: p.Trp691Cys, was found in family 2 and led to a Trp to Cys substitution. These two variants were absent in other healthy family members but were detected in the probands and affected relatives.

Two variants, c.1801G > A: p.Gly601Arg in LRP5 and c.633 T > A: p.Tyr211* in TSPAN12, were identified in families 3 and 4, respectively, and have been previously reported. The c.1801G > A missense variant in LRP5 resulted in the substitution of Gly with Arg at position 601, while the c.633 T > A missense variant in TSPAN12 led to a nonsense mutation at Tyr position 211, resulting in protein truncation.

Pathogenicity analysis

While sequencing validation was not conducted in the control population, the two variants in the LRP5 gene (c.4289delC: p.Pro1431Argfs*8 and c.2073G > T: p.Trp691Cys), as well as one variant in the TSPAN12 gene (c.633 T > A: p.Tyr211), were not identified in public databases such as dbSNP, ExAC, or 1000 Genomes. In contrast, the minimum allele frequency (MAF) of the c.1801G > A: p.Gly601Arg variant was 0.0000083.

In silico analysis methods SIFT, Polyphen-2, REVEL (scored 0.994, and 0.965 respectively), MutationTaster, and GERP +  + (scored 4.11, and 4.13 respectively) predicted two of the variants (c.2073G > T: p.Trp691Cys and c.1801G > A: p.Gly601Arg) to be disease-causing or damaging. The remaining two variants (c.4289delC: p.Pro1431Argfs*8 and c.633 T > A: p.Tyr211*) were predicted to be disease-causing or damaging by MutationTaster and GERP +  + (scored 4.53 and 5.68 respectively) (Table 2). All four variants exhibited co-segregation among affected family members. Notably, the newly identified missense variant c.2073G > T: p.Trp691Cys is evolutionarily conserved across a broad spectrum of species, extending from humans to zebrafish (Fig. 3B and C). Therefore, all these variations were considered pathogenic according to the ACMG standard.

Protein structure evaluation

The Ramachandran plot analysis of the three structures in Fig. 4 reveals that amino acid residues falling within the favored regions make up more than 90% of the entire protein. This suggests that the models' conformations conform to the fundamental principles of stereochemistry [52].

Fig. 4

The protein structures were modeled using AlphaFold2 and their conformations were evaluated via Ramachandran plots. A For the wildtype protein, 91.57% of amino acid residues were found in the favored regions. B In the p.Trp691Cys mutant, 90.08% of residues were in the favored regions. C The p.Pro1431Argfs*8 mutant showed 93.87% of its residues in the favored areas

Molecular dynamic simulation and data analysis

RMSD and RMSF have commonly used measures to assess the spatial variations of proteins in MD simulations. RMSD quantifies the overall difference of a molecule compared to a reference conformation. The RMSD values of the wildtype and two LRP5 mutations were illustrated in Fig. 5A. During the simulation, the p.Trp691Cys variant showed smaller conformational changes than the wildtype (P = 0.000, < 0.05, 95%CI) (Fig. 5E), while the p.Pro1431Argfs*8 variant did not show statistically significant differences in conformational changes compared to the wildtype (P > 0.05, 95% CI). In the end, the backbone atoms of all three proteins remained around 1.5 nm.

Fig. 5

RMSD, RMSF, Rg, HBNUM values comparisons between wildtype and two mutant proteins in MD simulation. A-D Evolution over time to show RMSD, RMSF, Rg, and HBNUM values of the wildtype (black) and the two mutant proteins (p.Trp691Cys in red, p.Pro1431Argfs *8 in blue). E Statistical analysis comparison of RMSD and RG values as bar graphs, depicting values with p < 0.05 (95% CI)

The RMSF analysis revealed distinct flexibility profiles among the wildtype and mutant proteins at specific residues (Fig. 5B) The wildtype's most flexible regions were Trp10-Ala27, Asp1318-Ile1387, Met1513-Asn1520, and Pro1606-Ser1615. In contrast, p.Trp691Cys displayed limited flexibility at Leu14-15, Asp1331-1333, and Ser1486-Thr1489; p.Pro1431Argfs8 had flexibility at Met1-Pro28 and Pro1261-Ser1385. Average RMSF values for the active sites fluctuated between 1.03–1.38 nm, 1.02–1.08 nm, and 1.42–1.53 nm for the wildtype, p.Trp691Cys, and p.Pro1431Argfs8, respectively. Maximum RMSF values were 1.87 nm, 1.16 nm, and 2.68 nm, respectively, implying reduced conformational changes in p.Trp691Cys and heightened local activity in p.Pro1431Argfs8, especially at Met1-Pro28 and Pro1261-Ser1385 (Fig. 6A). These findings suggest that the pathogenic variants may influence the flexibility of the protein backbone, which could affect the protein's structure and function. Especially evident in the videos of Additional file 1, the conformational changes of the three proteins over time were well-illustrated. Although the C-terminal region of the p.Pro1431Argfs*8 mutant exhibited heightened activity, it too reached a state of convergence by the end of the simulation.

Fig. 6

Comparison of active sites and H-bonds between the wildtype and mutant proteins. A. Changes in active sites in wildtype and two mutant proteins during simulation were marked in red. B. The snapshots extracted from the last frame of the MD simulation showed the intramolecular H-bonds changes in wildtype and mutants. The target residues along with H-bonds were marked in bright green

Rg values, which measure the distance or fluctuation from C-α atoms to the center of mass, were 4.34 nm, 4.25 nm, and 4.45 nm for the wildtype and two mutants, respectively (Fig. 5C). The C-α atoms in the p.Pro1431Argfs*8 mutant showed greater fluctuation (P = 0.000, < 0.05, 95%CI)) compared to the wildtype, while the p.Trp691Cys mutant exhibited less fluctuation (P = 0.000, < 0.05, 95% CI) (Fig. 5E). These results suggest that the p.Trp691Cys variant may increase the rigidity of the protein, whereas the p.Pro1431Argfs*8 variant may enhance its flexibility at the local conformation, especially in the regions of Met1-Pro28 and Pro1261-Ser1385.

HBNUM stands for the number of hydrogen bonds (H-bond), which are crucial interactions for stabilizing the spatial arrangement of proteins. It is evident from the results that the number of H-bonds in p.Trp691Cys is similar to that of the wildtype, while in p.Pro1431Argfs*8, which is a truncated protein due to the absence of C-terminus structures, the overall number of H-bonds is reduced compared to the wildtype (Fig. 5D). In terms of intramolecular interactions through H-bonds, as shown in Fig. 6B, Trp691 formed H-bonds with Ser700 in the wildtype, while the Cys691 of p.Trp691Cys formed H-bonds with Ala677, Ala679, and Ser700. In the p.Pro1431Argfs*8 mutant, the Pro1431 was transformed into Arg, which formed H-bonds with Tyr719, Glu721, and Thr737. The changes in H-bond distribution certainly contributed to the structural and conformational alterations.

To perform DSSP analysis, the last 50 ps of the simulation were selected, and the final frames of the MD simulations for the three proteins were obtained and compared with their respective initial structures. Additionally, a comparison was made between the 3D structures of the wildtype and the two mutant proteins in the final frame. All three ensembles were observed to converge to more stable conformations relative to their initial states, as depicted in Fig. 7A. In the final frame, both mutant proteins exhibited greater structural convergence and rigidity compared to the wildtype protein (Fig. 7C).

Fig. 7

Molecular Trajectories and Structural Comparisons of Wildtype and Mutant Proteins in MD Simulation. A Snapshots of molecular trajectories for the three proteins at both the initial and final frames of the simulation are superimposed and color-coded, with the initial frame marked in green and the final frame in red. B Quantitative analysis of eight principal secondary structures (coil, β-sheet, β-bridge, bend, turn, α-helix, 5-helix, and 3-helix) in the wildtype and the two mutants, p.Trp691Cys and p.Pro1431Argfs*8. C The figure displays snapshots of molecular trajectories at the final frame of the simulation, comparing the wildtype and two mutant proteins. The wildtype protein is marked in grey, while the two mutant proteins are denoted in green

Figure 7B presents the DSSP analysis results. Post-simulation, the wildtype protein showed a decrease in coil, β-sheet, and turn structures, while witnessing an increase in bend and 3-helix formations. The α-helix configuration remained constant, and a β-bridge structure emerged. Importantly, the 5-helix structure was absent throughout the entire simulation period.

Before the simulation, p.Trp691Cys had fewer β-sheet and turn structures than the wildtype but more bend, 3-helix, and α-helix structures. Post-simulation, turn and α-helix structures increased, surpassing the wildtype, while β-sheet remained stable. The β-bridge disappeared, and no 5-helix was observed at any stage.

For p.Pro1431Argfs*8, pre-simulation levels of β-sheet were similar to the wildtype, but turn structures were less abundant. Bend, 3-helix, and α-helix were more prevalent. Post-simulation, turn structures increased to wildtype levels; bend and α-helix decreased while 3-helix rose, still lagging behind the wildtype. A 5-helix structure appeared post-simulation. Throughout, the coil structure was consistently less abundant than the wildtype. The DSSP analysis is visually represented in Additional File 2.

Molecular interactions

In the docking simulation, the binding energies for wildtype-DKK1, p.Trp691Cys-DKK1, and p.Pro1431Argfs*8-DKK1 were -90.37 ± 2.35 kcal/mol, -135.78 ± 2.75 kcal/mol, and -118.69 ± 5.38 kcal/mol, respectively. Residue decomposition showed that hydrophobic and electrostatic interactions mainly drove binding, followed by nonpolar solvation-free energy. Both mutant proteins displayed larger negative values for van der Waals and electrostatic energies compared to the wildtype, indicating a stronger affinity for DKK1. Solvation factors for the mutants also suggested a more favorable binding environment with DKK1 (Table 3). These results indicate that both p.Trp691Cys-DKK1 and p.Pro1431Argfs*8-DKK1 have higher binding affinity compared to wildtype-DKK1 protein. Therefore, the pathogenic variants enhance the binding propensity of the protein with DKK1.

Table 3 Binding free energies and energy components predicted by MM/GBSA (kcal/mol)

Figure 8 illustrates the H-bond interactions between DKK1-C (residues 178–266) and the conserved third and fourth YWTD-EGF-β-propeller domains (residues 631–1246) of LRP5. Specifically, nine residues on DKK1 form H-bonds with nine residues on the wildtype protein.

Fig. 8

The 3D binding modes of the wildtype, p.Trp691Cys, and p.Pro1431Argfs*8 systems with DKK1 protein are shown, with a focus on the docking site. The LRP5 protein is represented in purple, while the DKK1 protein is depicted in a yellow–brown color

In the p.Trp691Cys-DKK1 complex, 11 residues from each protein are involved in H-bond formation. Similarly, in the p.Pro1431Argfs*8-DKK1 complex, 10 residues from each protein contribute to hydrogen bonding. These interactions provide insight into the differential affinities and stabilizing factors in the binding of wildtype and mutant LRP5 to DKK1.

Compared to the wildtype, both the p.Trp691Cys and p.Pro1431Argfs*8 mutants exhibit an increased number of H-bonds in their protein binding complexes, along with a change in the residue types that interact with DKK1 protein. Details on the specific H-bond interactions and bond lengths can be found in Additional File 3.

Alanine mutagenesis analysis

By comparing the mutated residues with alanine residues, the impact of point mutations on protein structure is elucidated. As shown in Additional File 4, the alanine variants at positions 691 and 1431 lead to a decrease in binding free energies and an increase in binding affinity compared to the wildtype. The binding free energy for the p.Trp691Ala is -134.75 ± 3.96 kcal/mol, while the binding free energy for the p.Pro1431Ala is -95.54 ± 3.18 kcal/mol. In contrast, the wildtype binding free energy is -90.37 ± 2.35 kcal/mol.