A list of 189 phosphatases, their phylogenetic classification, predicted substrate type and other descriptors were taken from (Chen et al., 2017), see (Supplementary Table S1). This included catalytic subunits from the following families: AP, alkaline phosphatase (4 genes), DSP, dual specificity phosphatase (40 genes) Myotubularin (15 genes), PTEN, phosphatase and tensin homolog (8 genes), PTP, protein tyrosine phosphatase (37 genes), Sac, Sac domain-containing phosphatase (5 genes), CDC25, cell division cycle 25 (3 genes), EYA, eyes absent family (4 genes), FCP, RNA polymerase II carboxy-terminal domain phosphatase (8 genes), NagD, NagD family of HAD-fold phosphatases (5 genes), HP1, histidine phosphatase 1 (12 genes) HP2, histidine phosphatase 2 (8 genes), PPM, PPM/PP2C family (20 genes), PPP (13 genes) and PPP-like phosphatases (2 genes), RTR1 (1 gene), other (3 genes).

Data was obtained via the secure Genomics England (GEL) research environment following information governance training as a members of the Genomics England Clinical Interpretation Partnership (GECIP): Enhanced Interpretation Domain, with approved project ID: 720—Towards understanding the aetiology of rare genetic disorders associated with altered phosphatase function. This provided access to data for 34,082 probands from 35,002 families (October 2022) analysed using the Illumina Starling pipeline and passing quality control parameters (Smedley et al., 2021). Patients were recruited to the study based on a residual unmet diagnostic need with a presumed underlying rare genetic disorder following standard clinical assessment and genetic investigations. These rare conditions were enriched in neurodevelopmental disorders, but also included disorders of growth, metabolism, the renal tract, skeletal and cardiovascular systems, vision, dysmorphic and congenital anomaly syndromes, tumour syndromes, and haematological and immunological disorders. After variant calling, filters were routinely applied to remove common variants, variants with no predicted impact and those that did not segregate with disease. Phenotypic descriptors, as assigned by the referring clinician, were recorded using Human Phenotype Ontology (HPO) terms. Variant details, associated patient phenotypic information and Genomic Medicine Centre (GMC) outcomes were extracted using the R LabKey package.

We used MARRVEL (Model organism Aggregated Resources for Rare Variant ExpLoration) via http://marrvel.org/human/batch/genes, to extract data from a range of human databases and assist in variant prioritisation (Wang et al., 2017). Additional information on existing syndromes was obtained from OMIM (https://omim.org/) and via literature searches. LOEUF (loss-of-function observed/expected upper-bound fraction) metrics were also compared across other predictive measures of intolerance to loss-of-function obtained from DECIPHER (https://www.deciphergenomics.org/). The predicted effect of genomic variants on known transcripts was determined using the Ensembl Variant Effect Predictor http://www.ensembl.org/Homo_sapiens/Tools/VEP (McLaren et al., 2010), with liftover of genome coordinates and annotation from GRCh37 to GRCh38 genome assembly using the UCSC LiftOver tool http://genome.ucsc.edu/cgi-bin/hgLiftOver, where necessary. Estimates of variant frequency in the general population were obtained from the genome aggregation database - gnomAD v3.1.2 or v2.2.1 (https://gnomad.broadinstitute.org/). To assist in the validation of candidate novel disease genes, additional patient cases for selected genes were identified from DECIPHER (Database of Genomic Variation and Phenotype in Humans Using Ensembl Resources) (https://www.deciphergenomics.org/) (Foreman et al., 2022). We curated high-level HPO terms for 100 kGP patients to enable comparison of patient phenotypes for selected variants across 100 kGP and DECIPHER. For example, “intellectual disability” was upscaled to “abnormality of the nervous system” (e.g. see https://hpo.jax.org/app/browse/term/HP:0000707).

To assess potential impact of variants identified for candidate disease genes we analysed Combined Annotation-Dependent Depletion (CADD) scores, which integrate diverse annotations for each variant, including allelic diversity, pathogenicity and annotations of functionality, into a single measure (Kircher et al., 2014; Rentzsch et al., 2021; Rentzsch et al., 2019). CADD scores were obtained by searching precomputed PHRED-like (-10*log10 (rank/total)) scaled C-scores for all possible human single-nucleotide variants (8.6 × 10⁹) via https://cadd.gs.washington.edu/snv.

Given the broad role of phosphatases in development and the availability of genomic data for undiagnosed rare disease patients, we considered it timely to systematically assess frequency of genetic variants in phosphatase genes and their possible association with rare-developmental disorders. We focused our attention on de novo changes identified in a cohort of patients recruited to the 100 kGP. We studied 189 genes encoding the human phosphatome (Figure 1), as defined by a recent genome-wide, cross-species survey of phosphatase catalytic subunit genes, which can be classified according to their structural folds and membership of subfamilies based on evolutionary origins (Chen et al., 2017). Amongst patients enrolled into the 100,000 Genomes Project (as of October 2022) we determined there were 220 rare de novo variants in 92 phosphatase genes (49% of the phosphatome), including 39% novel variants that are not present in the Ensembl database (see Supplementary Tables S1 and S2). The majority of MANE-select variants (Morales et al., 2022) were determined to be missense mutations (74%), with splice region variants (8%), premature stop (6%) and frameshift mutations (6%) being prevalent among the other classes (Figure 2A; Supplementary Table S2).

22 out of 37 phosphatase genes known to be associated with a rare disorder (as recorded in OMIM as of October 2022) were represented by at least one variant (60%). The most highly represented genes of this type were: PTPN11 (24 variants); PPP1CB (8 variants), PTEN (6 variants). Although most patients carrying variants in these genes had received a diagnosis, one pathogenic variant in each of PTPN11, PPP1CB and PTEN had not been reported as a diagnostic outcome. This is most likely because the recorded phenotypes did not prompt analysis of these genes, reinforcing the need for reinterrogation of the available genomic data in undiagnosed 100 kGP patients.

Although many of the variants we identified mapped to known disease genes, 61% of the total number of variants resided in genes that had not been previously associated with rare conditions, revealing a large potential for novel disease gene discovery (Figure 2B). The coverage of variants across different subfamilies of phosphatases varied considerably. Among multigene phosphatase groupings (i.e., not including single genes, such as Paladin), the maximum coverage of variants was in the PTP family (∼75% of genes/family with variants), the minimum coverage was in the CC2 enzymes (no variants), Figure 2C and Supplementary Table S1.

In total, we detected de novo variants in 70 genes not previously associated with a rare developmental disorder. To identify novel disease gene contenders we considered focusing on genes predicted to be intolerant to gene inactivation. This principle can be used to prioritise genes that have fewer deleterious loss-of-function (LoF) variants in the general population than would be expected by chance (Seaby et al., 2022) because natural selection removes variants that reduce survival and ability to reproduce (Karczewski et al., 2020). One metric to describe this is the loss-of-function observed/expected (O/E) upper-bound fraction (LOEUF). Genes in the first decile of LOEUF, with a score < 0.2, are amongst the most enriched for OMIM haploinsufficient disease genes (Karczewski et al., 2020). However, a survey of autosomal dominant inherited phosphatase disease genes in our dataset reveals that only three out of 13 (23%) disease genes (PTPN11, PPP1CB and PPP3CA) would meet this LOEUF criterion, with seven genes scoring 0.26–0.69 and three genes (PPM1D, PTPN22, PTPRQ) obtaining LOEUF scores of >1.0. Therefore, we adopted a filtering strategy that focused primarily on identifying genes with variants associated with matching phenotypes in unrelated patients.

The classification criteria and number of genes assigned to each of the categories are shown in Table 1. Class 1 genes included matches to patients in recent published datasets (after the end of 100 kGP patient recruitment, post 2018). Class 2 and 3 genes were identified as having variants in two or more unrelated kindreds in GEL with the same phenotype. We distinguished Class 2 and 3 genes using the upper bound of the O/E confidence interval LOEUF < 0.35 as a threshold. Variants were not considered in our analyses if an alternative diagnosis had been proposed by clinicians as recorded in the GEL research environment (GMC outcomes). This approach identified two class 1 genes: PPP2CA and PTPN4, for which multiple variants have recently been reported as being linked to neurodevelopmental disorders (Reynhout et al., 2019; Chmielewska et al., 2021), validating our approach. Pathogenic mutations in PPP2CA have been reported across the coding region of the gene (Reynhout et al., 2019). We uncovered missense mutations in exon two and three that add to spectrum of likely deleterious PPP2CA mutations associated with neurodevelopmental disorder and intellectual disability.

Having established that our approach was capable of identifying bone fide rare-disease genes, we focussed our attention on novel genes. After excluding variants in known disease genes, analysis of 134 variants using our selection criteria (Table 1) uncovered nine Class 2 and 3 genes (PTPRD, PTPRG, PTPRT, PTPRU, PTPRZ1, MTMR3, GAK, TPTE2, PTPN18) representing candidate rare developmental disorder genes (Table 2). Analysis of family structures for affected probands confirmed that all cases contained full parent offspring trios or families with more than three participants. Diagnostic yield amongst such families is approximately twice that for singletons for whom no other family member was recruited for whole genome sequencing (Smedley et al., 2021). To prioritise genes for further analysis we analysed variant abundance in the general population and potential pathogenicity. 67% of variants were not present in the general population database, consistent with them being deleterious changes. The highest variant frequency was 1.6 × 10⁻⁴. We assessed potential pathogenicity of variants using CADD scores, which rank variants relative to all possible substitutions of the human genome (8.6 × 10⁹). The median score for 100 kGP variants we identified (Table 1) is 25.1, placing them in the top 1% most deleterious possible substitutions.

To provide further evidence that gene perturbation is linked to the associated phenotypes, we interrogated cases reported in the Decipher database (Wright et al., 2015) to find additional de novo variants in the selected genes with matching high level HPO terms. Seven out of nine genes had variants associated with matching phenotypes in Decipher. Although there is some potential for patient overlap in Decipher and 100 kGP, all but one variant was distinct between the two databases and therefore come from unrelated patients. To rank the candidate genes for future follow up, we used a points-based system to assign a confidence score in the candidate genes (Table 2). Together, our findings suggest that phosphatases represent promising candidate genes for a range of undiagnosed developmental disorders and in particular, neurodevelopmental syndromes, which are relatively underexplored.

We have focused here on phosphatase catalytic subunits. However, it is important to note that the majority of PPP catalytic subunits complex with regulatory subunits that direct these enzymes’ specific roles. Consequently, non-catalytic subunits may also represent a source of promising disease gene candidates. Indeed, nine genes encoding regulatory subunits of PP1 and PP2A are known to be associated with heritable disorders (Vaneynde et al., 2022). However, systematic analysis of the non-catalytic subunits poses a challenge because of the large number of disparate proteins that interact with PPP catalytic subunits. This is particularly true of PP1 which is thought to bind more than 180 PP1-interacting proteins (Verbinnen et al., 2017), with diverse functions ranging from molecular chaperones to signalling molecules and transmembrane receptors and channels. Moreover, because many PPP binding proteins act as multivalent binding or scaffold proteins, non-catalytic subunits often have additional functions unrelated to their phosphatase-binding properties. For example, SARA, Spinophilin and various A-kinase anchoring proteins (AKAPs), bind to a variety of PPP catalytic subunits and exert diverse regulatory functions by controlling the assembly, trafficking and localisation of signalling complexes containing multiple components (Bennett and Alphey, 2002; Wong and Scott, 2004; Shaw and Filbert, 2009; Rozes-Salvador et al., 2018). Additional roles for PPP regulatory proteins that are not exclusively dependent on the phosphatase catalytic activity are continually being uncovered, highlighting the complexity of their functions (Sakaguchi et al., 2022). These features of the non-catalytic subunits complicate the ability to predict phosphatase dependent and independent effects of variants in these genes. In future, careful scrutiny of variants in genes encoding non-catalytic subunits, accompanied by analysis of potential phenotypic overlap with the relevant catalytic subunits will begin to address this issue. The findings that we report here sets the foundation for ongoing efforts in this direction.

Foremost of the novel disease gene candidates we have identified is PTPRD, for which there are several lines of evidence supporting its role in neurodevelopmental disorders. In situ hybridization analysis has revealed that murine PTPRD is expressed in the specialized regions of the brain including the hippocampal CA2–CA3 region, thalamic reticular nucleus, piriform cortex, olfactory bulb, olivary nucleus and spinal motor neurons (Mizuno et al., 1993; Sommer et al., 1997; Schaapveld et al., 1998; Shishikura et al., 2016). Moreover, studies into the fundamental cell and developmental roles of PTPRD in mice have shown roles for PTPRD in neural precursor and cortical development (Tomita et al., 2020), axonal growth and pathfinding (Uetani et al., 2006). Clinically, genome wide association studies have linked PTPRD variants with several neural disorders, including: 1) rare copy number variation with autism spectrum disorders (ASDs) (Pinto et al., 2010), attention-deficit/hyperactivity disorder (ADHD) (Elia et al., 2010), and obsessive-compulsive disorder (OCD) (Gazzellone et al., 2016); 2) Single nucleotide polymorphisms with schizophrenia (Li et al., 2018); and, 3) non-coding variants associated with restless leg syndrome (Schormair et al., 2008; Yang et al., 2011). Our identification of de novo PTPRD protein coding changes in patients affected by motor defects and speech and learning difficulties raise the prospect of there being a neurodevelopmental syndrome associated with altered PTPRD function. Notably, Ptprd knockout mice also show defects in motor development as well as learning and memory (Uetani et al., 2000). This and other models may therefore offer an ideal opportunity to assess functional effects of the variants we have identified.

PTPRT is also a prominent neurodevelopmental disease gene candidate. Expression of PTPRT RNA and protein is enriched in the brain (McAndrew et al., 1998; Uhlen et al., 2015), and in early development PTPRT transcripts are distributed throughout the murine brain and spinal cord (McAndrew et al., 1998), suggesting a role in CNS function. Functional studies in mice indicate that these roles are likely to include hippocampal neurogenesis, synaptogenesis and dendritic arborization of hippocampal neurons (Lim et al., 2009; Lee, 2015; Lim et al., 2020). Notably, an individual with intellectual disability and a de novo PTPRT variant has previously been reported (Schuurs-Hoeijmakers et al., 2013). Further cases will help validate the role of PTPRT in neurodevelopmental disease. We are now building case series and undertaking functional analyses to confirm genotype-phenotype relationships for PTPRD, PTPRT and other genes we have identified.

Our gene discovery approach relies heavily on finding a match between phenotypes of unrelated patients with de novo variants in a given gene of interest. However, since the cohort is enriched with similar phenotypes there is a risk of false positives arising. Conversely, due to the variability and non-standardised application of HPO terms used to describe similar phenotypes, there may also be false negatives. Our previous analysis across two genomic medicine centres (North West Coast GMC and Greater Manchester GMC, n = 3,212 index cases) suggests that diagnostic yield is poor in cases with limited or non-specific phenotypic information, reflected in a low recorded number (≤ 15) of HPO terms (E Lyulcheva-Bennett 2020; personal communication). The occurrence of inaccurate or incomplete phenotype records highlights the need to collaborate with clinicians to obtain comprehensive and up-to-date phenotypic data.

The majority of novel variants in candidate disease genes that we discovered are not found in the general population database, gnomAD. Those variants that were found at low frequency might be non-pathogenic, or pathogenic with incomplete penetrance, modified by cis-regulatory variation, or associated with adult-onset disease (Gudmundsson et al., 2021). We are unable to distinguish between these possibilities at the current time. We used CADD scores to obtain a preliminary measure of potential pathogenicity. Almost all identified genes have at least one variant with a CADD score >25, which is comparable to well characterised pathogenic lesions in known phosphatase disease genes (e.g., PTPN11-T73I, CADD 25.6; PPP1CB-P49R, CADD 26.7).

Haploinsufficiency metrics have gained popularity in the field of disease gene discovery as they can be used to focus the scope of genomic data reanalysis and to inform disease etiology. We did not use these metrics for candidate disease gene identification because known phosphatase disease genes span a broad spectrum of haploinsufficiency tolerance. In addition, natural selection is blind to phenotypes that do not affect the ability to reproduce. Therefore, genes associated with later-onset phenotypes, or phenotypes not impacting on reproduction, exhibit much weaker intolerance to inactivation and will most likely not be picked up when applying haploinsufficiency constraint-based metrics. Moreover, the relationship between genetic constraint and molecular pathogenesis is likely to be complex. While the pathogenic mechanism of highly constrained genes might be assumed to be loss-of-function, this is not always the case. For instance, pathogenic gain-of-function mutations in the constrained PTPN11 gene (LOEUF 0.14) are a known cause of Noonan’s Syndrome. On the other hand, loss-of-function PTPN11 mutations give rise to an allelic disorder (Noonan’s Syndrome with Multiple Lentigines), with overlapping clinical features (Solman et al., 2022). Similarly, variants in genes that are not predicted to be haploinsufficient, can nevertheless be pathogenic, for instance by acting in a dominant-negative fashion, as is the case for C-terminally truncating mutations in PPM1D (LOEUF 1.1) (Jansen et al., 2017).

With the aspiration to sequence five million more genomes in the United Kingdom as part of routine clinical practice, and with increasing numbers of patients undergoing whole genome sequencing globally, we anticipate that reiterative gene-centred analysis of abundant genomic and phenotypic data will help to identify further phosphatase gene variants for functional validation. This should go hand-in-hand with clinical gene discovery clinics re-examining unsolved cases. Our current focus is on expanding the scope of our analysis to include de novo variants affecting the phosphatase regulatory subunits and uncover novel recessive phosphatase disorders. The reanalysis of genomic data provides the first step towards systematically delineating the aetiology of rare genetic disorders associated with altered phosphatase function. Many of the known developmental disease genes were first identified in model species based on their role in fundamental developmental processes. Consequently, experimental model systems provide a rich resource with which to dissect mode of action, establish disease causality, and determine variant pathogenicity. Improved information sharing among clinicians, genomic scientists and biologists will assist in efforts to aid phosphatase disease gene discovery and delineate genotype-phenotype relationships.