The Division of Medical Genetics, Fondazione IRCCS Casa Sollievo della Sofferenza, includes a regional referral center for hereditary bone metabolism disorders. The laboratory receives requests for genetic analysis from different Italian tertiary care centers. For the current study, cases were recruited from five pediatric (RF, TF, CG, ES, PV), three endocrinology (ASS, AG, RMM/SC), and four bone metabolism (FP/CB/AS, FC, SC/GG, MCe) referral centers for adults, two medical genetic (EG/SM, MC) units. From March 2017 to August 2022, 30 index individuals were enrolled for the study after preliminary biochemistry assessment confirmed the ALP value lower than the minimum range and measured at least twice and the VitB6 level higher than the maximum range (7–18 μg/l). Exclusion criteria were the assumption of drugs affecting bone turnover (i.e., bisphosphonates, denosumab) and/or diseases associated with low ALP activity/Vit B6 value (i.e., hypomagnesaemia, hypozincemia, celiac disease, Wilson disease, hypothyroidism, hypoparathyroidism, assumption of multivitamin supplements).

Three additional patients (#21, #22, and #23) underwent a next-generation sequencing (NGS) multigene panel analysis sequencing for hereditary disorders leading to early-onset osteoporosis/bone fragility and were included in the final survey after the identification of an ALPL variant: (i) subject #21 was a male newborn with fractures in the absence of trauma at 30 days of life with a suspicion of “shaken baby” syndrome; (ii) subject #22 was a boy of 2 years and 10 months, with an initial diagnosis of EDS/bone fragility; (iii) patient #23, was an adult man of 65 years, with persistent osteopenia and fractures. Unexpectedly, subjects (#21) and (#23) showed a normal (repeatedly) ALP value ( Table 1 ).

Totally, 48.4% of the subjects were from South Italy, 27.2% from the Islands (Sardinia and Sicily), and 24.4% from North-Middle Italy. The sex ratio was F/M = 1.72 and median age = 51.8 +/- 19.4 years. Clinical, biochemical, and anamnestic data were collected.

DNA was extracted from blood withdrawal, and the classic Sanger method was performed on 25 (75.7%) out of the 33 patients: DNA was PCR amplified, and amplicons were purified (ExoSAP-IT, Affymetrix, Thermo Fisher) and sequenced (Big Dye Terminator Cycle Sequencing Kit v 1.1, ABI3130XL Sequencer, Thermo Fisher) for the 12 coding exons ( Supplementary Material 1 for primers and PCR cycling conditions).

Eight patients (24.3%) were screened through a SureSelect gene panel (Agilent Technologies, USA) designed to capture known genes associated with calcium metabolism disorders and bone fragility, including ANO5, ALPL, B3GALT6, B4GALT7, BMP1, COL1A1, COL1A2, CREB3L1, CRTAP, DKK1, EN1, FAM46A, FKBP10, FKBP14, IFITM5, LRP5, MBTPS2, P3H1, P4HB, PLOD1, PLOD2, PLOD3, PLS3, PPIB, SEC24D, SERPINF1, SERPINH1, SLC39A13, SPARC, SP7, TMEM38B, WNT1, and YY1AP1.

Targeted fragments were sequenced on a MiSeq platform (Illumina, USA) with MiSeq Reagent kit V3 300-cycle flow cells. Data analysis was performed considering the frequency, impact on the encoded protein, conservation, and expression of variants using distinct tools [ANNOVAR (9), dbSNP (10), 1000 Genomes (11), ExAC (12)]. Selected variants were interpreted according to the American College of Medical Genetics and Genomics/Association for Molecular Pathology (ACMGG/AMP) (13). Variants were confirmed by Sanger sequencing and segregation investigation performed on available relatives.

MLPA Kit (P484-ALPL, MRC Holland) was initially applied on negative subjects and then enlarged to the whole cohort, and it was carried out following the manufacturer’s instructions.

The level of conservation of the protein sequence was evaluated for the novel missense variants. From the Ensembl Genome Browser (14), the TNSALP protein orthologous sequences of the following organisms were downloaded: Homo, Gorilla, Pan, Mus, Bos, Xaenopus, Tetraodon, Danio. Then, Clustalw (15) was used to create the multi-alignment.

Atomic coordinates of TNSALP protein were obtained through AlphaFold v.2.0 (16). We retrieved the ALPL protein sequence in FASTA format from the Ensembl Genome Browser (14) and submitted it to ColabFold (17), which could predict the complex folding using the specialized AlphaFold-multimer model.

Then, molecular dynamics (MD) simulation was employed to refine the obtained models. Both the monomer and TNSALP dimer were embedded into a simulation box, extending up to 12 Å, using the tleap module of AmberTools21 (18), and solvated with the TIP3P water model and an appropriate number of Na+ and Cl- counter ions. The Amber ff14SB force field was applied to amino acids. Each model was subjected to energy minimization using steepest descent, followed by conjugate gradient methods. Then, they were gradually heated and equilibrated for approximately 5 ns, by time steps of 1 fs. A 10-Å cutoff was used for non-bonded short-range interactions, whereas long-range electrostatics were treated with the particle-mesh Ewald method. Temperature and pressure were maintained at 300 K and 101.3 kPa, respectively, using the Langevin dynamics and Piston method (19).

Following these preliminary steps, we implemented the Gaussian-accelerated molecular dynamics (GaMD), an enhanced molecular dynamic sampling method that works by providing a harmonic boost potential, which follows a Gaussian distribution, to smooth the system potential energy surface, thus causing the decrease of local energy barriers and accelerating the transitions between low-energy states. The boost potential can be applied in a dual-boost scheme: (i) a dihedral potential boost and (ii) a total potential boost. Maximum, minimum, average, and standard deviation values of the system potential were obtained from an initial ~2 ns conventional simulation with no boost potential, followed by a further GaMD equilibration of ~50 ns in which the boost potential was updated every 1.6 ns, thus reaching the equilibrium. Thus, each system was simulated for 100 ns using a timestep of 2 fs.

Finally, the stability of multiple TNSALP variants, p.(Ser181Leu), p.(Tyr28Asp), p.(Val95Met), p.(His267Pro), p.(His472Arg), p.(Met219Ile), p.(Ala179Thr), p.(Ala33Val), p.(Arg223Gln), (p.Arg391His), p.(Asp109Glu), p.(Gln207Pro), p.(Glu235Gly), p.(Glu298Lys), p.(Glu84Lys), and p.(Glu88Lys), and the double mutant p.(Arg136His)-p.(Phe290Leu), was investigated thermodynamically through the BuildModel and Stability function implemented in the FoldX algorithm (20). It was run with standard parameters. The same workflow was employed to evaluate the stop loss p.(*525Argfs*11) but using ColabFold to model the structure of both the monomer and dimer containing an additional residue fragment.

The original values of ALP were measured for non-uniform ranges; therefore, to make them comparable, the “ALP percentage deviation” = [(lower ALP value (for the corresponding reference interval)) - ALP)/lower ALP value (for the corresponding reference interval) x 100] was calculated. Analysis was performed with R (packages: stats, ggplot2, plyr, dbplyr, dplyr) (21–25).

Apart from early-onset cases or the ones that were resolved through the NGS-differential diagnosis (#21, #22 and #23), the remaining patients were firstly admitted for clinical problems related to bone pain, myalgia, fractures, nephrolithiasis, or early-onset or post-menopausal osteoporosis screening. This biased sample selection influenced the overall clinical spectrum of the cohort; instead, around 48% of patients claimed osteoporosis/low bone mass whereas 27% showed fractures at different sites. However, it is important to highlight that the whole cohort was then enrolled only after the biochemistry showed low ALP and high VitB6 ( Table 2 ).

From a total of 30 index individuals, we identified one or two single-nucleotide variants in 19 and a single-exon deletion (exon 2, Figure 1 ) in an additional individual, for a total of 20 cases (66.7%). To this cohort, we added three additional cases (#21-23), who received the diagnosis of HPP after the NGS multigene panel testing identified a likely pathogenic ALPL variant. Totally, 21 cases were heterozygous and two resultant compounds were heterozygous, a finding compatible with autosomal recessive inheritance (cases #5 and #10). Of the 25 identified mutated alleles (i.e., 21 heterozygous plus 2 compound heterozygous), we identified 22 different variants (21 single-nucleotide variants and 1, exon 2, deletion). Accordingly, 18 were private whereas 3 (i.e., c.262G>A, c.327C>A, and c.668G>A) recurred in 2 apparently unrelated pedigrees/index individuals. The 21 single occurrences were distributed as follows: 16 were missense, 3 frameshift, with the presumed insertion of a premature stop codon, 1 stop loss with the presumed insertion of 11 additional residues, 1 stop gain. Eight variants out of 21 (38%) were novel since they were not reported in the official database (5) (last release, March 2023). The Franklin online genomic variant interpretation tool (26) was used to determine the pathogenicity level of the eight novel changes, which resulted in one pathogenic, five likely pathogenic, and two VUS, following the ACMGG guidelines (13) ( Table 1 ).

There were 12 out of 23 mutated cases (52%) found to be familial, although for others, a clinical in-depth study or segregation analysis was not available, or yet concluded at the time of this report ( Figure 2 ). Interestingly, 10 (out 30, 33.3%) subjects resulted to be negative, although at least for one case (#29) familiarity for low ALP value and clinical symptoms was reported.

Due to the different assays and normal ranges used in the tertiary cares, we decided to plot the percentage deviation of the ALP value from the minimum value versus the presence/absence of the ALPL mutation. The correlation between the presence of a deleterious ALPL variant and the corresponding ALP serum value was, then, resolved in the boxplot in Figure 3 : mutated patients showed an ALP value that was from 11 up to 49% lower than the minimum.

Other biochemical parameters (Ca, PTH, P, VitD, and Cr) were collected, but no association was found with the presence/absence of an ALPL variant, nor with the ALP value.

Several SNPs were found among positive and negative patients ( Supplementary Material 2 ). p.Arg152His is defined in ClinVar (8) as benign, with a MAF = 1.5%, in gnomAD. However, in our cohort, it was found in three subjects (9%), two of them carrying a pathogenic variant, with an overall frequency ~6 times higher than expected ( Supplementary Material 2 ).

The multi-alignment of the Alpl ortholog clearly showed that the novel missense variants affect residues conserved up to the amphibians, or up to the fishes, strongly supporting that they could have a deleterious effect on the protein function ( Figure 4 ).

FoldX computed the total energy of ALPL dimers, wild-type and mutated, as a proxy of their overall stability and the van der Waals inter-residue clashes, as energy penalization factors. The free energy calculations are reported in Table 3 and Figure 5 . The difference in free-energy ΔΔG (ΔGmt – ΔGwt) for the majority of the ALPL amino acid changes was in the range to classify them as destabilizing whether only one or both monomers in the dimeric complex were mutated. In detail, p.(Asp109Glu), p.(His472Arg), p.(Met219Ile), p.(Ser181Leu), p.(Ala179Thr), p.(Arg391His), p.(Glu84Lys), and p.(Tyr28Asp) were observed as destabilizing whereas p.(Gln207Pro), p.(His267Pro), p.(Ala33Val), and p.(Arg223Gln) and the double mutant p.(Arg136His)-p.(Phe290Leu) could be classified as highly destabilizing. Conversely, p.(Glu88Lys) could be classified as stabilizing.

p.(Glu298Lys) exhibited stabilizing effects when only one monomer was mutated; however, this alteration was compensated in the double-mutated dimer; on the contrary, p.(Val95Met) showed a destabilizing effect in the mutant monomer but a stabilizing effect with two mutations. Finally, p.(Glu235Gly) did not exhibit any structural impact on the ALPL protein.

A different approach was applied to the stop loss p.(*525Argfs*11), where the additional residue fragment in the protein structure could have a putative impact on the overall folding. Comparing the stability of both the monomer and the dimer with the wild type, we observed how the mutated protein was highly unstable in its monomeric form (ΔΔG = 555.75 kcal/mol), whereas the dimeric configuration was slightly more stable than the wild-type dimer (ΔΔG = -126.52 kcal/mol). It has to be noted that three mutants (p. (Tyr28Asp), p.(Ala179Thr), and p.(Glu235Gly)) were not part of the mutant collection found in our cohort, but they were previously studied (27, 28) and here included as internal controls.

In our cohort of suspected 30 HPP cases, we observed a mutation rate of 66.7% (20/30), sensibly higher if compared with previous studies (43%, 47%) (32, 33). Such a difference might lay on heterogeneity of patients’ selection, technical variability, and genetic background. Instead, Tenorio et al. (32) recruited 83 subjects on the basis of their clinical manifestation and low ALP activity, but the PLP value was not considered. Conversely, the cohort of 72 patients reported by Jandl et al. (33) was selected by low ALP activity, higher PLP, and clinical manifestation but consisted of only adult subjects (41.3–68.5 years). In the first case, it can be argued that adding a higher PLP level as an additional inclusion criterion could improve the case selection and the mutation rate. To confirm this assumption, we used, as internal control, a different cohort of subjects with low ALP activity but with normal or (presumably normal) not provided VitB6 value. Screening of the ALPL gene on these subjects was negative (data not shown), indicating that although hypophosphatasemia (with normal VitB6 serum level) is a not infrequent biochemical issue, it would not be related to the TNSALP gene. Concerning the Jandl et al.’s study, it cannot be excluded that the older age of the sample size may have influenced the mutation rate, perhaps due to the overlap with general osteoporosis.

In our survey, 10 patients with low ALP and high VitB6 values, one also having familial hypophosphatasemia, resulted to be negative. With regard to the clinical picture of these patients, whether the phenotype could be softened or worsened by the presence of alleged pathogenic variants lying in genetic regions not covered by the screening (such as deep intronic regions or 5′ and 3′UTR) or compensatory mechanisms or the action of genetic modifiers is far from being elucidated. In this regard, we kindly point out that our NGS panel contained the whole genomic (intronic and regulatory regions) of the ALPL gene, in order to identify deep intronic splicing variants, which, however, had never been found.

Interestingly, patients’ enrollment was driven by biochemistry values, regardless of the clinical spectrum that, as reported in Table 2 , was highly variable, ranging from low bone mass/bone pain or just myalgia and actually reflecting the paradigmatic different HPP clinical variability. Thus, our work can confirm the utility of biomarkers in case of blurred clinical presentation, since only the ALP and VitB6 values were confirmed as the most suggestive clues with a strong correlation (up to 65% of cases) with the presence of a genetic lesion ( Figure 3 ). In our previous work, we also suggested the predictive value of PEA (34); however, this metabolite, often along with VitB6, is not included in the routinely biochemistry panel, and, then, no PEA data were available at the time of this report.

Literature data recorded the p.(Gln207Pro) variant once in a severe autosomal recessive perinatal form in compound heterozygosity with p.(Arg71Pro) in a baby girl who died 3 h after birth (2). Patient #21 carried the same 207P variant and was affected by an autosomal dominant perinatal (30 days of life) form featured by severe fractures. However, the mother and grandmother, both carriers of the same variant, had an asymptomatic form with borderline ALPL levels. The p.(Glu88Lys) variant was found once in compound heterozygosity with the p.(Pro37Leu) in an 11-month-old boy affected by nephrocalcinosis and alterations of the fontanelle (35). Both our patients #7 and #23 carried the 88Lys in heterozygosity and #23 had nephrocalcinosis, since he was 10 years old. Whether the 88Lys variant is specifically accountable of the nephrocalcinosis requires further studies. It should be noted that cases #21, #22, and #23 were not initially included in the survey, as ALP values for #21 and #23 were normal. However, the use of NGS in the differential diagnosis helped to clarify, in both, an uncovered HPP form.

Interestingly, 2 subjects, #6 and #9 carrying the same variant, the p.D109E, both came from Sardinia and, although they denied any apparent kinship, a possible ancestor effect cannot be excluded.

About the eight novel amino acid changes, the ALPL database has two different records for residue Glu84 (#18): the first was 84Val found in a family with autosomal dominant odonto-hypophosphatasia also characterized by rickets-like signs (36). The second, 84Asp, was found as a compound heterozygote in a 6-year-old girl with a childhood form. Despite no functional test being done on these variants, the authors suggested that both variants at Glu84 would bear a DNE, since the residue is located at the homodimer interface (36): this assumption was confirmed by our bioinformatic investigation.

No data can be inferred from the literature about the effect of the 267Pro (#15) or 472Arg (#2) mutants. Conversely, with regard to the Phe290 variant, Mornet et al. proved that this residue, together with Glu235, Glu291, and Asp306, coordinates the binding of the calcium atoms (37): it can therefore be deduced that the 290Leu (#10) variant could severely impair the enzymatic activity. Our modeling study confirms that all three mutants destabilize the homodimer structure.

About the stop gain p.(Tyr85*) (#19) and both the frameshift mutations, p.(Ser310Profs*28) (#17) and p.(Arg391Profs*14) (#11), that, presumably, introduce a premature stop codon, the mutated proteins are likely to undergo proteasomal degradation (38), or, if translated, they are totally inactive, due the loss of (more than) 50% of the canonical protein structure. In contrast, p.(*525Argextfs*11) (#13) does not interrupt the coding frame, as it leads to correct translation of the whole protein with presumed 11 additional residues. Although in the absence of functional data, the effect of this supplemental tail could not be considered detrimental for the enzymatic activity; instead, the mutant does not seem to bear the DNE since the phenotype of patient #13 (and of her brother, Figure 2 ) overlaps with classical osteoporosis with a borderline-low ALP activity. In our computational modeling, the mutant homodimer p.(*525Argextfs*11) appears more stable than the monomer, so we cannot exclude that the destabilized structure is partially mitigated. This is in contrast with data reporting that a similar mutation, p.(Leu520Argfs*86), that adds more than 80 residues at the C-terminal tail was proven to lack GPI binding with a consequent severe reduced enzymatic activity (39).

There were 22 common intronic (nr 14) and exonic (nr 8) variants identified: 17 have a MAF in gnomAD > 0.05, whereas the remaining were rarer ( Supplementary Material 3 ). No correlation was attempted with the presence of fractures, nor with biochemical values, as previously suggested (40, 41), since the distribution of the SNPs in the whole cohort (positive, negative, low, or normal ALP patients) was not different (data not shown).

p.(Arg152His) was found in three patients (out of 33, 9.3%) with a frequency six times higher than that reported in gnomAD (MAF = 1.5%): in two cases (#7 and #19), it was identified together with a pathogenic ALPL variant. Whether 152His could actually play the role of the phenotype modifier or it is in linkage disequilibrium with an actual pathogenic variant located in genetic regions not analyzed would be a hypothesis to shed light on. Moreover, it cannot be excluded that its frequency is higher in the Italian population, although there is no evidence in public databases.