Digenic inheritance of subclinical variants in Noonan Syndrome patients: an alternative pathogenic model?

Noonan syndrome (NS) is an autosomal-dominant disorder with variable expressivity and locus heterogeneity. Despite several RAS pathway genes were implicated in NS, 20-30% of patients remain without molecular diagnosis, suggesting the involvement of further genes or multiple mechanisms. Eight patients out of 60, negative for conventional NS mutation analysis, with heterogeneous NS phenotype were investigated by means of target resequencing of 26 RAS/MAPK pathway genes. A trio was further characterized by means of whole-exome sequencing. Protein modeling and in silico prediction of protein stability allowed to identify possible pathogenic RAS pathway variants in four NS patients. A new c.355T>C variant in LZTR1 was found in patient 43. Two patients co-inherited variants in LRP1 and LZTR1 (patient 53), or LRP1 and SOS1 genes (patient 67). The forth patient (56) carried a compound heterozygote of RASAL3 gene variants and also an A2ML1 variant. While these subclinical variants are singularly present in healthy parents, they co-segregate in patients, suggesting their addictive effect and supporting a digenic inheritance, as alternative model to a more common monogenic transmission. The ERK1/2 and SAPK/JNK activation state, assessed on immortalized lymphocytes from patients 53 and 67 showed highest phosphorylation levels compared to their asymptomatic parents. These findings together with the lack of their co-occurrence in the 1000Genomes database strengthen the hypothesis of digenic inheritance in a subset of NS patients. This study suggests caution in the exclusion of subclinical variants that might play a pathogenic role providing new insights for alternative hereditary mechanisms.

Introduction

Noonan syndrome (NS) is a RASopathy with autosomal-dominant inheritance and prevalence of 1:1000–2500 [1]. NS is characterized by variable expressivity of the phenotype, including craniofacial features such as hypertelorism, downslanting of the palpebral fissures, ptosis, low-set, posteriorly rotated ears, and webbed or short neck [1-3]. More than 80% of individuals with NS exhibit cardiovascular involvement, most frequently congenital heart diseases, pulmonary valve stenosis and hypertrophic cardiomyopathy [4]. Other clinical manifestations include cryptorchidism, bleeding disorders, mild neurocognitive delay and pectus deformity, and an increased risk of developing myeloproliferative disorders [3, 5]. Given the considerable variable expressivity, some NS adults are diagnosed after the birth of an affected child [6, 7].

NS shows significant locus and allelic heterogeneity, since variants causing disease, affecting gene function, involve different genes of the RAS/MAPK (mitogen-activated protein kinase) signal transduction pathway: PTPN11 (in 40–50% of the patients), SOS1 (10–20%), RAF1 (3–17%), RIT1 (9%), and the less common KRAS, NRAS, BRAF, SHOC2, MAP2K1, CBL, LZTR1, SOS2, RRAS, CDC42 and A2ML1 [8, 9] (1–5%) [9]. The A2ML1 variants have been reported to be associated to few NS patients [10] even if other authors interpreted one of the two previously reported A2ML1 pathogenic variants, of unknown significance (VUS) [11]. Differently, a significant number of data evidenced the relevance of functional role of LZTR1 in RAS signaling considering the effects of its functional variant in NS patients [9, 12], Interestingly, the finding of biallelic LZTR1 pathogenic variants reported in 12 families with NS children, supports an autosomal recessive inheritance besides the more frequent dominant inheritance pattern with different implications for the NS pathogenesis [13]. Moreover, a few NS and NS–like cases have been reported to show Copy Number Variations (CNV) encompassing NS-associated genes [12, 14]. However, the pathogenesis remains unknown in about 20–30% of the patients, pinpointing the need for identifying new genes or mechanisms responsible for NS pathogenesis [5, 15]. Interestingly, we have previously reported a double heterozygous patient carrying both an inherited SOS1 and a de novo RAF1 variants [7]. As the family carriers of that SOS1 variant alone exhibited a subclinical phenotype, we have hypothesized that the co-expression of additional subclinical mutant effectors of the RAS pathway may have an addictive effect thus contributing to the pathogenesis. By means of NGS target approach, we performed the resequencing of RAS/MAPK pathway genes on eight patients with a clinical diagnosis of NS and negative for the conventional NS mutation analysis. In three unrelated patients, we have identified novel affecting function variants in the recently described LZTR1 NS-gene and in two new identified genes, LRP1 and RASAL3, never associated with NS. Furthermore, the co-occurrence of sub-clinical variants in two RAS-pathway genes in the two probands prompted us to propose digenic inheritance as an alternative pathogenic mechanism occurring in some NS patients, negative to conventional genetic analysis.

Materials and methods

Patients and cell samples

Informed consents for genetic studies and publication of case reports were obtained from all enrolled individuals or parents. Blood samples of these patients had been collected in the laboratories of the University of Milan mainly from two hospitals, Fondazione Policlinico San Matteo in Pavia and Istituto Auxologico Italiano in Milan The treatment of human subjects complies with the Declaration of Helsinki, and the Research Ethics Committee of IRCCS Istituto Auxologico Italiano (ID number 08C607_2006) approved the study. Epstein Barr Virus (EBV) immortalized PBMCs, were provided by Galliera Genetic Bank from patients’ lymphocytes [12] (Supplementary File).

Sequencing analysis and quantitative real-time PCR (qPCR)

Constitutional DNA and total RNA were extracted from whole blood according to standard procedures. PCR was carried out following standard procedures using specific oligonucleotides (Table S1). The PCR products were sequenced and resolved on an automated ABI-3130xl DNA genetic analyzer (Thermo Fisher).

Total RNA (500 ng) was reverse transcribed using the iScriptTM cDNA Synthesis Kit (Bio-Rad Laboratories Inc., Barkeley, CA, USA) and specific oligonucleotides allowed to amplify A2ML1 gene (Table S1), TBP gene was used as a housekeeping control (Supplementary File).

NGS analysis and protein modeling

The enrichment of regions of interest was performed using Agilent SureSelect XT technology. Target resequencing was performed on eight patients’ samples by Illumina MiSeq platform. Subsequently, whole-exome sequencing (WES) was performed on a trio, as the proband was found uncharacterized in the previous analyses. For each sample, a MiSeq run was performed. Raw reads of NGS data are available in NCBI Short-read Archive (SRA, https://www.ncbi.nlm.nih.gov/sra/) under accession number PRJNA607236.

The detailed bioinformatics analyses applied to NGS data are reported in Supplementary File. To achieve multiple sequence alignment, necessary for protein modeling, protein sequences were obtained from Uniprot and the CLUSTAL O servers were employed. For protein modeling, multiple protein sequences were retrieved from Uniprot and protein models were obtained through the Swissmodel Target-Template Alignment module. The final alignments were carried out by means of ClustalO application (Supplementary File).

Array comparative genomic hybridization (aCGH) analysis

The aCGH experiments were performed using a commercial SurePrint G3 ISCA CGH + SNP Microarray Kit, 4 × 180 K slides (Agilent Technologies Inc., Santa Clara, CA, USA) in accordance with the manufacturer’s instructions.

Immunoblotting

Proteins were extracted in 7 M urea, 2 M thiourea, 4%CHAPS, 30 mM Tris, 1 mM PMSF (Sigma Aldrich), and 10 µl/ml phosphatase inhibitor cocktails A and B (Santa Cruz Biotechnology). Blots were incubated with anti-SAPK/JNK, anti-phospho SAPK/JNK (Cell Signaling Technology), anti-ERK1/2 and anti-phosphoERK1/2 (ThermoFisher Scientific) rabbit polyclonal antibodies (Supplementary File).

Statistical analysis

The ratio of phosphorylated over the total SAPK/JNK and ERK1/2 band intensity after western blot was evaluated by applying the ANOVA with post-hoc Tukey HSD (Honestly Significant Difference) tests (n = 3, p value < 0.01 or p value < 0.05). The co-occurrence probability of the two SOS1 and LRP1 variants was calculated considering their minor allele frequency (MAF) in both the whole 1000genomes dataset and in the European subgroup (503 individuals). We also carried out the Fisher exact test (http://vassarstats.net/odds2x2.html), and the Peto Odd Ratio (OR) was calculated for the events of co-occurrence of the two variants in both our cohort of patients and in the European subgroup of 1000genomes.

Results

Patients

Eight patients (Pts 43, 49, 53, 55, 56, 62, 67, 69) out of a cohort of 60, negative after conventional sanger sequence analysis, with or suggestive of NS phenotype, have been molecularly characterized. The clinical features of the above patients are shown in Table 1, with the exception of Pt 62, who has been meanwhile diagnosed to have mosaicism for a small supernumerary marker chromosome derived from chromosome 19 and described by Recalcati et al. [16] as patient 1. Whenever possible, for each of the clinical feature reported in Table 1 the corresponding item in the HPO (Human Phenotype Ontology) phenotype vocabulary [17, 18] has been indicated.

Table 1

Clinical featuresa of NS patients .

Patient n., Gender	43, F	49, M	53, M	55, F	56, M	67, F	69, F
Growth parameters at	10 m	31 gw	40 gw	38 gw	35 gw	37 gw	na
OFC (percentile)	cm 45 (50th)	na	cm 32 (1st, −2.57 SD)	na	cm 35 (95th)	na	na
Length (percentile)	cm 66 (1th, −2.28 SD)	na	cm 43 (0th, −4.07 SD)	na	cm 44 (10th)	cm 49 (64th)	na
Weight (percentile)	g 7700 (10th)	g 1025 (3rd, SGA)	g 2540 (0th, −2.59 SD, SGA)	g 2600 (8th, SGA)	g 3117 (92nd, LGA)	g 2400 (9th, SGA)	na
Growth parameters at	11 y 3 m	3 y	21 y 7 m	9 y 7 m	7 y 9 m	9 y 1 m	2 y 5 m
OFC (percentile)	cm 55 (75th)	Absence of macrocephaly	cm 55 (9th)	Microcephaly	cm 49.5 (0.4th)	cm 50 (0.4th)	cm 46.5 (0.4th)
Height (percentile)	cm 143.8 (>25th)	cm 94.6 (25th)	cm 145.5 (−4.98 SD)	cm 119 (−2.69 SD)	cm 120.5 (10–25th)	cm 94 (−7 SD)	cm 77.6 (−3.51 SD)
Weight (percentile)	Kg 37 (<25th)	Kg 12.9 (10th)	Kg 48 (−3.11 SD)	Kg 22 (−2.19 SD)	Kg 20.85 (5–10th)	Kg 9 (−9 SD)	Kg 9.1 (−3.59 SD)
Development	Normal	Normal	delayed (HP:0001263)	na	Delayed (HP:0001263)	Severely delayed (HP:0012736)	Normal
♦ID (HP:0001249) (degree)	−	/	+ (moderate) (HP:0002342)	+	+ (moderate with severe language delay) (HP:0002342, HP:0006863)	+ (severe with psychomotor regression between 6–9 y) (HP:0010864, HP:0002376)	/
Craniofacial
♦Low posterior/trident hairline (HP:0002162)	−	−	+	−	+	−	−
♦Hypertelorism (HP:0000316)	−	−	+	+	−	−	+
♦Epicanthic folds (HP:0000286)	+	−	+	+	+	+	+
♦Downslanting palpebral fissures (HP:0000494)	+	+	+	−	+	−	−
♦Ptosis of eyelids (HP:0000508)/thick eyelids (HP:0030939)	+	−	+	+	+	−	−
♦Posteriorly rotated ears (HP:0000358) with thick elix (HP:0000391)	+	−	−	−	+	−	+
♦Myopathic face (HP:0002058)	+	−	−	−	+	−	−
♦Short/broad neck (HP:0000470; HP:0000475)	−	−	+	−	+	−	−
♦Webbed neck (HP:0000465)	−	−	−	−	+	−	−
Ocular
Strabismus (HP:0000486)	−	−	−	+	+	−	−
Skeletal
♦Superior pectus carinatum (HP:0000917) and inferior pectus excavatum (HP:0000915)	−	+	−	+	+	−	−
♦Broad thorax (HP:0100625) with wide-spaced nipples (HP:0006610)	−	−	+	−	+	−	+
Cubitus valgus (HP:0002967)	−	−	+	−	−	−	+
Clinobrachydactyly (HP:0040019; HP:0001156)	−	−	+	−	−	−	−
CHD (HP:0001627)
♦Pulmonary valve stenosis/branch pulmonary artery stenosis (HP:0001642)	+	−	−	−	−	−	+
♦Atrial septal defect (HP:0001631)	−	−	−	−	+	+	−
♦Ventricular septal defect (HP:0001629)	−	+	−	−	−	+	−
♦Other CHD	Slight asymmetric hypertrophy of inteventricular sept (HP:0005144)	−	−	Persistence of fetal left superior vena cava (HP:0005301)(diastolic heart murmur)	Aortic arch hypoplasia	−	−
♦Cryptorchidism (HP:0000028)	/	−	+	/	+	/	/
Lymphatic dysplasias
♦Chylous effusions of the pleural space and/or peritoneum (HP:0002202, HP:0012281)	−	−	−	−	+	−	−
♦Cystic hygroma in early pregnancy (HP:0010878)	−	−	−	−	+	−	−
Epidermal
Cafe au lait spots (HP:0000957)	+	−	−	−	−	−	−
Nevi or lentigines (HP:0003764; HP:0001003)	+	−	−	−	−	−	−
Permanence of fetal finger (HP:0001212) and toe pads	+	−	−	−	−	−	−
Other	Narrow palate (HP:0000189)	Low-set ears (HP:0000369)	GH deficiency (HP:0000824); previous epilepsy (rolandic type, absences); thoraco-lumbar scoliosis (HP:0002944); generalized hirsutism (HP:0002230)	Sparse hair (HP:0008070), hypoplastic supra-orbital ridges (HP:0009891), blepharophimosis (HP:0000581), depressed nasal bridge (HP:0005280), low-set ears (HP:0000369), high palate (HP:0000218); long and tapering fingers (HP:0100807; HP:0001182), overlapping toes (HP:0001845); delayed tooth eruption (HP:0000684)	Polyhydramnios (HP:0001561), hydrops fetalis (HP:0001789); sensorineural deafness (HP:0000407); slight lumbar scoliosis (HP:0004626); back hirsutism (HP:0009889)	Generalized hypotonia (HP:0001290) and hypotrophy (HP:0003202), but normal muscular strength; GH deficiency (HP:0000824); onset (7 y) of episodes of hyperkinetic movements (HP:0002487) at upper limbs, likely not epileptic; diagnosis of painful joints and panhypopituitarism (HP:0000871) at 9 and 11 y	Dysplastic and low-set ears (HP:0000377; HP:0000369), micrognathia (HP:0000347)
Genetic tests	46,XX (at CVS)	nd	46,XY; wild-type for array CGH as well as for FGFR1, FGFR2, FGFR3 mutational screening, and FRAXA	na	46,XY; wild-type for FRAXA; GRCh37: chr1:g.153207564_153598790del (inherited from the mother), chr7:g.29148038_29180543del (inherited from the father)b; normal muscle and fibroblast respiratory chain analyses; normal dosage of VLCFA, plasma and urinary amino acids, urinary organic acids, free and esterified carnitine;	46,XX; wild-type for array CGH as well as for LMNA, PTEN, GLUT1, MECP2, CDKL5, FOXG1, MEF2C, SCN1A, PCDH19, STXBP1, and UBE3A sequencing; normal dosage of muscle enzymes, plasma and urinary amino acids, urinary organic acids, plasma lactic acid, pyruvic acid, ammonia, and sialotransferrine electrophoresis	Wild-type for array CGH

♦features that define the diagnostic criteria of NS [3], + present, − absent, / not applicable, na not available, nd not done, gw gestational week, y years, m months, SD standard deviation, SGA and LGA small and large for gestational age, VLCFA very long chain fatty acids, chr chromosome, del deletion.

aFeatures frequently seen in NS but not exhibited by our patients were not reported in the table.

bThese Copy Number Variants (CNVs) have been interpreted as likely benign. Whenever possible, for each of the clinical feature within the brackets has been reported the corresponding item in the HPO (Human Phenotype Ontology) phenotype vocabulary.

After the NGS target resequencing, the clinical phenotypes of the patients were meticulously re-assessed by means of the same clinician. Patients 43, 53, 56 and 67 could be re-evaluated together with their family members, while for patients 49 and 55 only the clinical form filled in during the patients’ enrollment could be revised. Patients 43, 53 and 56 have a clinical diagnosis of NS based on the criteria defined by van der Burgt et al. [3] and exhibit the typical NS face dysmorphology. Similarly, Pt 69 could be suggestive for NS diagnosis. Patient 67 instead had been candidate for screening in the spectrum of RASopathies due to relative macrocephaly, hypotonia, poor growth and growth hormone deficiency, congenital heart defect, and after the exclusion of the other hypothesized clinical conditions. From their clinical records, it can be hypothesized that Pt. 55, exhibiting hypoplastic supra-orbital ridges and blepharophimosis, may have a different diagnosis from NS. Similarly, Pt 49 exhibited very few of the typical NS craniofacial features and a CHD not specific to NS (Table 1). All of the parents and siblings that could be reviewed (families of Pts 43, 53, 56 and 67, Fig. 1) were confirmed unaffected. Families’ history was unremarkable with the exception of that of Pt 56: he was the second born fourth child from Moroccan consanguineous parents, as his paternal grandfather, II-3, and the mother, I-1, of his maternal grandfather, II-4, were siblings (Fig. 1d). The first pregnancy (IV-3) of his mother was reported to interrupt at the seventh month for polyhydramnios; the second pregnancy (IV-4) was reported to end at the fourth month of gestation without an apparent reason. The relatives-clinical records were not available. Both parents were karyotyped to rule out balanced chromosomal rearrangements. Moreover, the paternal cousin IV-2 was reported to be affected by hypertrophic cardiomyopathy; the paternal cousin IV-1, who was born from consanguineous parents too (III-1 and III-2 are cousins), by slight intellectual disability. Unfortunately, their genotypes, although requested in Morocco, were not carried out.

Fig. 1

Pedigrees of NS patients and segregation analysis of the identified variants.

Filled symbols, subjects with NS/NS spectrum; light-dotted filled symbols, reported to be affected by slight Intellectual disability (d, IV-1) or hypertrophic cardiomyopathy (d, IV-2); asterisk, members phenotyped and genotyped in the study; Pt patient, SA spontaneous abortion.

NGS target resequencing

To identify the affecting function variants associated with the clinical condition of the enrolled patients, we designed a target sequencing panel including genes conventionally associated with RASopathies, genes recently associated with NS and genes belonging to the RAS/MAPK pathway and never associated with RASopathies, for a total of 26 genes (Table S2). We generated 2.9 million reads per sample on average. After duplicate removal, we obtained an average of 2.8 million reads mapped on NCBI human reference genome (build GRCh37). The mean depth was 747×(ranging from 384 to 899), with more than 99% of the targeted regions covered by NGS reads in each sample (Table S3). Only 317 variants passed all the filtering steps imposed by our pipeline (e.g., low-depth, strand-bias etc.) (Supplementary File). By means of Annovar, the annotation analysis revealed 312 single-nucleotide polymorphisms (SNPs) and five deletions. The majority of SNPs occur within intronic regions (114), while the remaining are located in the 3′UTR or downstream (99), in 5′UTR or upstream (18), in intergenic regions (20) and in exons (61). Among the coding variants, 16 and 45 are missense and synonymous respectively (Table S4). Aiming at selecting coding non-silent rare variants in 1000genomes database (1000g2015aug_eur), we applied a stringent filter, according to a MAF < 0,05, obtaining eight rare variants in five genes (Table 2). To reduce false-positive results, we applied a combinatorial approach by comparing different prediction tools for assessing the possible effects of nonsynonymous SNPs (snSNPs) (Table S5). All the selected variants were validated by Sanger sequencing. No affecting function variants were identified in patients 49, 55, 62 and 69. Chromosomal unbalances were ruled out by aCGH analysis. Proband 43 was found to be a carrier of the novel missense variant c.355T>C (p.(Tyr119His)) in the exon 4 of LZTR1 (Fig. 1a, Table 2). This de novo variant was predicted deleterious by the majority of in silico tools applied (Table S5). The novel variant c.1430C>T (p.(Ala477Val)) in the exon 13 of LZTR1, was identified in heterozygosis in Pt 53 (Table 2; II-2, Fig. 1b) and was paternally inherited. The variant was never reported (Table 2) and was predicted damaging by different in silico tools (Table S5). Moreover, three maternally inherited LRP1 missense variants were identified in heterozygosis: c.136C>T (p.(Arg46Trp)), c.5977C>T (p.(Arg1993Trp)), c.8591G>A (p.(Arg2864His)) (Table 2; II-2, Fig. 1b). This gene encodes a member of the low-density lipoprotein receptor family of proteins, and interacts with FGF receptors and CBL, mediating its ubiquitination activity involved in RAS pathway regulation. All these variants were reported as deleterious, in particular the c.5977C>T (p.(Arg1993Trp)) variant was predicted damaging by all considered tools (Table S5). The patient 67 showed the paternally inherited variant c.650C>T (p.(Ala217Val)), rs1800127 (Table 2; II-4, Fig. 1c) in the exon 6 of LRP1 predicted as tolerated (Table S5). This patient also showed the maternally inherited SOS1 variant c.1964C>T (p.(Pro655Leu)), rs56219475, that was also carried by the unaffected sister II-3 (Fig. 1c). The in silico predictions on the effect of this variant were discordant (Table S5), but the variant was previously reported as benign [6]. Given these results, we hypothesized an addictive effect of the two variants that co-segregate in the proband. At this purpose, we interrogated 2504 individuals included in the 1000genomes database (phase 3). The two LRP1 and SOS1 variants display a frequency of 0.00698 and 0.00319 in the 1000genomes database, while in our cohort we detected 0.0199 and 0.0119, respectively. Furthermore, we assessed if these two variants are never present in the same individual. We used Fisher exact test to compare the frequency of co-occurrence in the 1000 genomes EUR (N = 503 individuals, no occurrences) and in our sample set (N = 8 patients, 1 co-occurrence), and this difference is statistically significant (p value = 0.015, Peto OR = 5.50). Given the frequency of the two variants, the probability of co-occurrence in our small sample set (8 patients) is as low as p value = 0.00189, suggesting that this event should not be expected in the considered cohort, indicating that co-occurrence of LRP1 and SOS1 in our cohort with a frequency 0.125 is not a random event. Furthermore, to verify whether this event is due by chance, we evaluated in the 1000genomes database the co-occurrence of any two rare variants in all the genes included in our NGS panel. According to the previous analysis, we retained the nonsynonymous variants with a MAF < 0.05% and predicted to be damaging by at least two predictor tools. Following the identification of the individuals carrying more than two variants, we found that only 16 subjects, out of 2504 individuals (frequency 0.0064), shared the two variants SOS2 c.1448G>A (p.(Ser483Asn)) and A2ML1 c.1799T>C (p.(Val600Ala)). The novel missense variant c.2882A>G (p.(Asp961Gly)) in exon 24 of A2ML1 (Table 2) was found in heterozygosis in patient 56 (Table 2; IV-6, Fig. 1d). The variant was shared with both unaffected consanguineous parents (III-5 and III-6) and with the unaffected sister (IV-5). In addition, II-3 and II-4, who were reported to be healthy, are obligate carriers. All these findings suggest that c.2882A>G (p.(Asp961Gly)) is not sufficient to cause the phenotype (Fig. 1d). We, then, investigated whether differential allelic expression (DAE) in whole blood as well as in epithelial tissue could be responsible for the incomplete penetrance. No significant difference of A2ML1 specific alleles expression was observed comparing the proband to his sibling (data not shown).

Table 2

Variants identified by RAS panel target sequencing and WES.

Patient	Gene	Exon	dbSNP	Nucleotide variation	Protein variation
43	LZTR1 NM_006767.3	4	–	c.355T>C	p.(Tyr119His)
53	LZTR1 NM_006767.3	13	–	c.1430C>T	p.(Ala477Val)
	LRP1 NM_002332.2	2	rs764409398	c.136C>T	p.(Arg46Trp)
		37	rs141826184	c.5977C>T	p.(Arg1993Trp)
		54	rs143285614	c.8591G>A	p.(Arg2864His)
67	LRP1 NM_002332.2	6	rs1800127	c.650C>T	p.(Ala217Val)
	SOS1 NM_005633.3	12	rs56219475	c.1964C>T	p.(Pro655Leu)
56	A2ML1 NM_144670.5	24	–	c.2882A>G	p.(Asp961Gly)
	RASAL3 NM_022904.2	13	rs199734851	c.1963G>A	p.(Gly655Ser)
		7	rs58123634	c.751C>G	p.(Leu251Val)

For full details see Tables S4 and S5.

Whole-exome sequencing of trio 56

To elucidate the genetic bases of NS in family 56, we hypothesized the involvement of genes different from those belonging to the RAS pathway and we performed WES in the trio. At total of 49.11 million reads per sample was obtained; after duplicates removal, and an average of 47.03 million reads were mapped on the NCBI human reference genome (build GRCh37). On average, the mean depth was 37X with more than 96% of the targeted regions covered by NGS reads in each sample (Table S6). All the selected variants were validated by Sanger sequencing. The condition of compound heterozygosis was detected for two variants in RASAL3 gene, encoding a protein with pleckstrin homology (PH), C2, and Ras GTPase-activation protein (RasGAP) domains: the c.1963C>A (p.(Gly655Ser)) variant (rs199734851) in the exon 13 of RASAL3 was paternally inherited, while the c.751C>G (p.(Leu251Val)) (rs58123634) in exon 7 of the same gene was maternally inherited (Table 2; Fig. 1d). In silico predictions of the deleterious effect of these variants were reported in Table S5.

Protein modeling and in silico prediction of protein stability of deleterious nonsynonymous variants

Sequences of genes with substitutions in encoding regions were analyzed. In detail, the protein structure of A2ML1, LZTR1, RASAL3, SOS1 and LRP1 were modeled (see Supplementary File for details). The obtained model of A2ML1 allowed us to evaluate the c.2882A>G(p.(Asp961Gly)) identified variant. The residue substitution localizes in the protein core and is not directly involved in post-translation modification (PTMs) or disulfide bridges. The 961 residue seems to have a role in proteases-binding or in their recognition and the substitution of a charged residue with a small apolar one may modify the protease binding activity. This substitution may be destabilizing because it is very close to the typical thiol ester sequence (969–973) [19] involved in the thiol ester bond (970–973) (Fig. 2a). Interestingly, this sequence is distinctive of the α-Macroglobulin, is involved in protein reactivity [19] and it is localized on the opposite side of the bait domain on a large binding cavity; the bait domain is generally cleaved by proteases, in order to induce a conformational change and entrap the protease itself [19].

Fig. 2

Structural predictions of A2ML1, LZTR1, RASL3, and SOS1.

A2ML1 structural prediction, substituted residue, Asp961Gly, is in blue stick and functional regions are in evidence, thiol ester sequence in pink and bait region in red (a). LZTR1 model and focus on Tyr199His substitution, in gray wild type (wt) residue and in green the model of Tyr in position 119 (b). RASAL3 model and focus on Gly655Ser substitution, in gray wt residue and in white Ser in position 655 (c). SOS1 crystal structure and focus on Pro655Leu substitution, in gray wt residue and in light-blue Proline in position 655 (d).

The LZTR1 model allowed us the evaluation of the two identified variants (Table S7): the substitution c.355T>C (p.(Tyr119His)) (Fig. 2b) that is predicted to destabilize the protein structure while the substitution c.1430C>T (p.(Ala477Val)) is located outside the modeled protein (ending in position 325) portion and is predicted neutral while (Table S7). Noticeably, His119 in its neutral form can establish π-π stacking interactions with the neighbor His120, therefore possibly reducing the flexibility of His121, which interestingly is a catalytic residue.

Concerning to RASAL3, all the three stability prediction tools reported c.1963G>A(p.(Gly655Ser)) variant as destabilizing for the protein structure (Table S7). The substituted residue localizes on the protein surface (Fig. 2c), within a loop, at the end of the functional RAS-GAP domain and therefore the substitution could be involved in the binding interface with RAS destabilizing this docking.

Considering the variant identified in SOS1 (Table S7), we assessed that the substitution c.1964C>T (p.(Pro655Leu)) is localized on the protein surface, in a loop, within the RAS-GEF domain (Fig. 2d). Since Proline residues are highly functional for protein structure also in loop conformations, this variant may be destabilizing for the protein structure.

The new candidate NS LRP1 gene, it encodes a large protein of about 4500 residues, mainly composed by LDL-receptor repeated domains [20] and by EGF-like domains, not equally distributed along the sequence (Fig. 3e). Due to the protein length, we were not able to obtain a single model for the whole structure, but we obtained four models not super-posable including residues 1 to 3332 (Fig. 3 a–d) (Supplementary File for details). All of the four identified substitutions of LRP1- c.136C>T (p.(Arg46Trp)), c.650C>T (p.(Ala217Val)), c.5977C>T (p.(Arg1993Trp)) and c.8591G>A (p.(Arg2864His))—were predicted to be destabilizing for the protein structure as indicated by negative difference of free energies of wild type and mutant structures (Table S7) and were modeled. The first substitution, c.136C>T (p.(Arg46Trp)) (Fig. 3a), is localized in a partially folded region and this substitution may influence cysteine 65, is involved in a disulfide bridge with Cys48. The substitution c.650C>T (p.(Ala217Val)) (Fig. 3b) is conservative and localizes on the receptor surface, therefore the involvement in the NS pathogenesis cannot be ascribable to protein structure destabilization. The substitution c.5977C>T (p.(Arg1993Trp)) (Fig. 3c) displays an increased steric hindrance in a packed region at the internal side of the YWTD-beta propeller. Moreover, substituting a positively charged residue with an aromatic nonpolar aminoacid destroys the high number of hydrogen bonds and salt interactions, which involve Arginine. This condition induces higher flexibility of the involved protein portion. Furthermore, a larger sidechain may have effects on the glycosylation of residue Asp1995, necessary for protein structure. The variant c.8591G>A (p.(Arg2864His)) (Fig. 3d) is localized in a disordered region, but in the nearby of a YWTD-beta propeller. This type of substitution may be classified as not destabilizing because a positively charged residue is substituted with another positively charged residue, but with a lower isoelectric point. However, the variant may influence the disulfide bridge of residue 2865 with Cys 2884.

Fig. 3

Assembly of LRP1 models and focus on LRP1 variants.

Each color identifies a partial model, starting and ending residue of each model is labeled. In the focus the detail of Arg46Trp substitution, in red wild type (wt) Arg and in white Trp variant (a), of Ala217Val, in red wt Ala and in white Val variant (b), of Arg1993Trp, in orange wt Arg and in white Trp variant (c) and of Arg2864His, in blue wt Arg and in white His variant (d). The domains organization is displayed (e), green pentagon represents EGF domains, green rhombus represents LDL-receptor class A domain, while in pink are represented YWTD domains.

Extracellular signal-regulated kinase (ERK) and Stress-activation protein kinase/Jun-amino-terminal kinase (SAPK/JNK) in patients 53 and 67

We carried out specific assays to verify the activation of both ERK1/2 and SAPK/JNK, thanks to the availability of immortalized lymphoblastoid cell lines obtained from the PBMCs of both trios 53 and 67, showing a co-segregation of more hypomorphic variants. In fact, the evaluation of ERK1/2 phosphorylation increase is a common assay to validate a new NS variant, but more recently a SOS1 variant has been observed to activate SAPK/JNK instead of ERK in an NS patient [7]. While the levels of ERK1/2 phosphorylation, assessed by immunoblotting, were not significantly different in patients 53 and 67 compared to their healthy parents and unrelated controls (Fig. 4), SAPK/JNK was activated in both NS patients. The ratio of the phosphorylated over the total SAPK/JNK band intensities was calculated for each sample. As expected, the anti-SAPK/JNK antibody revealed two bands at ~46 (p46) and ~54 KDa (p54), respectively (Fig. 5a, b). Of them, only p54 resulted phosphorylated in both patients (Fig. 5d). A statistically significant increment (ANOVA + Tukey, n = 3, p value < 0.01) of phospho-SAPK/JNK was detected in sample 67 compared to 67M, 67P and unrelated control (Fig. 5e) and in sample 53 compared to 53M, 53P and unrelated control (Fig. 5f). In addition, for sample 53, a second band, p54B, was detected both in total and in phosphorylated SAPK/JNK blots (Fig. 5b).

Fig. 4

Activation ERK1/2 assay.

Representative immunoblot images of stress-activated protein kinase (ERK1/2) (a), and phosphorylated ERK1/2 (b). Histograms show the ratio of phosphorylated over the total ERK1 (p44) band intensity (O.D. optical density, mean ± SD) in patient 53 and 67 compared to their parents and unrelated controls 524 and 548 (c). Histograms show the ratio of phosphorylated over the total ERK2 (p42) band intensity (O.D. optical density, mean ± SD) in patient 53 and 67 compared to their parents and unrelated controls 524 and 548 (d). Phosphorylation levels were not significantly different in patients 53 and 67 compared to their healthy parents and controls, applying the ANOVA + Tukey HSD tests (n = 3, significance threshold: p value < 0.05).

Fig. 5

Activation SAPK/JNK assay.

Representative immunoblot images of stress-activated protein kinase/Jun-amino-terminal kinase (SAPK/JNK) (a, b), and phosphorylated SAPK/JNK (c, d). Histograms show the ratio of phosphorylated over the total SAPK/JNK band intensity (O.D. optical density, mean ± SD) in patient 67 (e) and 53 (f) compared to their parents and unrelated controls 524 (e) and 548 (f). Significant differences are computed applying the ANOVA + Tukey HSD tests (n = 3, significance threshold: plain lines = p value < 0.01, dotted lines = p value < 0.05).

Discussion

Molecular diagnosis remains unknown at least in 20–30% of NS patients [5], therefore new genes and genetic mechanisms are expected to play a role in the disease. NGS approach led to the identification of new NS/NS-like genes [4, 5, 21], all encoding members of RAS/MAPK pathway. Most of them have been identified in less than 10 individuals [4], suggesting that all RAS/MAPK genes should be investigated as potentially implicated in the pathogenesis of NS, and more in general of RASopathies. Accordingly, we designed an extensive NGS target resequencing panel. A novel de novo missense variant c.355T>C (p.(Tyr119His)) was identified in LZTR1 in proband 43, who exhibits a typical NS phenotype with mild severity, since her development, intellectual skills and growth were unremarkable (Table 1). LZTR1 encodes a protein member of the BTB-kelch superfamily. Recently, a significant amount of data referring to functional role of LZTR1 in RAS signaling and the impact of these variants on protein function, have been reported [8, 22]. Yamamoto et al. identified rare variants of LZTR1 among which (p.(Tyr119Cys)) was reported as a de novo event occurring with NS phenotype [23]. As the substitution here described is not conservative, involves the same aminoacidic residue previously described [23] and may affect a nearby catalytic residue, we propose it as a new pathogenic NS-associated variant. Another novel variant in LZTR1 (NM_006767.3): c.1430C>T, predicted to be deleterious, but inherited from an unaffected father, was identified in proband 53. This variant could not be modeled and was not associated with NS phenotype if inherited alone. Nonetheless, in proband 53, we observed that the co-occurrence of LZTR1 (NM_006767.3): c.1430C>T variants with the three LRP1 variants NM_002332.2: c.[136C>T;5977C>T;8591G>A], all inherited from the unaffected mother (Fig. 1b), segregates with NS phenotype. This complex condition suggests a digenic inheritance model. LRP1 is involved in intracellular signaling functions as well as in lipid metabolism. It is a physiological modulator of the platelet-derived growth factor (PDGF) signaling pathway and can directly interact with the ubiquitin E3-ligase CBL, already associated with NS [24, 25]. This evidence suggests a Sprouty-like role for LRP1 in modulating ubiquitination and endocytosis of the tyrosine-kinase (Trk) receptors such as PDGFRbeta [26, 27]. Moreover, LRP1 regulates the A2ML1 activity by promoting its internalization, thus maintaining its clearance [28]. Interestingly, also the proband 67 paternally inherited the LRP1 (NM_002332.2): c.650C>T variant, which was not associated with NS phenotype if inherited alone. In addition, the proband 67 also maternally inherited the SOS1 (NM_005633.3): c.1964C>T missense variant (Fig. 1c). This variant is not rare in the population and has been previously reported as benign [6], however the Proline residue is highly functional for the protein structure, therefore the substitution could be destabilizing for the structure. All the identified LRP1 variants do not seem to be sufficient to determine the manifestation of NS phenotype, as we observed in both cases the co-occurrence with LZTR1 (Pt. 53) and SOS1 variants (Pt. 67). While the LRP1 and SOS1 variants detected in patient 67 are singularly reported in genomes from 1000genomes data set, interestingly, their co-presence does not occur. Both the calculation of co-occurrence probability and the application of Fisher exact test and OR, indicate that the co-presence of the identified subclinical LRP1 and SOS1 variants is expected as a very rare event. The identification of this condition in one patient of a cohort of eight individuals (frequency 0.125), together with the very low frequency of the co-occurrence of two any rare variants in the 1000 genomes population (frequency 0.0064), is consistent with its pathogenic significance, invoking the hypothesis of a digenic model in NS. Digenic inheritance of subclinical variants could partially account for pathogenic causes in 30% of diagnosed NS patients, explaining the lack of identification of affecting function variants in this NS subgroup [5]. This hypothesis is consistent with the increase of phospho-SAPK/JNK in the PBMCs cell lines of both patients in comparison to those of their parents, who carried only one of the two mutated genes present in the probands. We observed no significant changes in ERK1/2 phosphorylation, which is the functional marker commonly considered to validate NS variants [5]. Nevertheless, Longoni and colleagues reported on the increasing of phospho-SAPK/JNK after overexpression of SOS1 R497Q variant detected in the father and grandfather of a double heterozygous NS patient, carrying both SOS1 and RAF1 variants [7], where both proband’s relatives showed subclinical NS traits. It is, therefore, possible that SAPK/JNK pathway, which shares some common effectors with ERK1/2 pathway, might also be involved in NS pathogenesis. We speculate that the co-occurrence of variants in both LRP1 and in SOS1 or LZTR1 may alter RAS pathway upstream causing hyperactivation, likely by an addictive effect, of either ERK1/2 or SAPK/JNK in a specific tissue fashion [13].

Patient 56 shared the A2ML1 (NM_144670.5): c.2882A>G variant with all the unaffected family members. Segregation of the variant in the family suggests that it may be not associated with NS, nevertheless, our model predicted a possible effect on A2ML1 protein function. Indeed, family member IV-1 (Fig. 1d), who is affected by slight intellectual disability, may be homozygous for the above variant in A2ML1. Because A2ML1 is a protease inhibitor and binds to LRP1, regulating activation of the MAPK/ ERK cascade, as in the case of LRP1, we hypothesized that also this variant in A2ML1 could contribute to alter the physiological modulation of RAS/MAPK pathway and display NS insurgence in association with variants occurring in other genes. At this purpose we performed WES analysis in the trio and found that Pt.56 was a compound heterozygote for RASAL3 (NM_022904.2): c.[1963C>A]; and c.[751C>G] variants, which were maternally and paternally inherited, respectively. The finding, as well as a low probability of loss-of-function intolerance (pLI = 0.00) [29], is suggestive for RASAL3 being involved in NS pathogenesis according to an autosomal recessive mechanism. Indeed, the experience of polyhydramnios in the third trimester during the first pregnancy (IV-3, Fig. 1d) of Pt.56’s mother is highly suggestive of NS and of the likely existence of compound heterozygosity for RASAL3 pathogenic variants at least in IV-3 aborted fetus. Noteworthy, Pt.56, who exhibited the typical NS gestalt and had a history of polyhydramnios and hydrops fetalis during pregnancy, showed developmental delay, severe speech and language delay associated with sensorineural deafness, and moderate intellectual disability (Table 1). RASAL3 is a Ras GTPase-activating protein, i.e., RasGAP [30, 31] still little characterized. For instance, NKT cells from liver of RASAL3-deficient mouse treated with α-GalCer have been shown to increase Erk phosphorylation, suggesting that RASAL3 negatively regulates Ras/Erk signaling [31]. Despite this evidence may be encouraging about our finding, we believe that it should have been fundamental to investigate how a predicted loss of function of a RasGAP predominantly expressed in cells of hematopoietic lineages [30, 31] can result in a severely/fully expressed NS RASopathy. Unfortunately, PBMCs samples of family 56 were not available, therefore it was not possible to investigate ERK1/2 or SAPK/JNK phosphorylation levels. Furthermore, RASAL3 may also perform a Ras-independent function in neuronal signal transduction. Indeed, despite in the literature RASAL3 expression is studied in the lymphoid cells, according to BrainSpan (https://www.brainspan.org/) and GTEx (https://www.gtexportal.org/home/) databases, the gene is brain expressed prenatally (average level of expression = 0.186) and in infancy (average level of expression = 0.208). Moreover, it keeps on being expressed, at a higher level, during adult life (average level of expression = 2375.5) [32]. Further functional studies will be needed to confirm our findings. Here we provide a first evidence that digenic inheritance may be an alternative mechanism involved in the etiopathogenesis of NS and RASopathies, as already hypothesized [7, 33]. Moreover, we have identified two novel genes, LRP1 and RASAL3, possibly involved in the etiopathogenesis of NS/RASopathies. In particular, RASAL3, whose affecting function variants have been identified in a family from a population with a high rate of inbreeding, represents, after LZTR1, the second gene associated with autosomal recessive inheritance in NS so far reported [34]. Here we have identified the presence of variants in major genes of the RAS/MAPK pathways in a subgroup of NS patients selected from a cohort 60 with heterogonous NS phenotypes, but negative for conventional NS mutation analysis. To gain more confidence in the proposed alternative hereditary mechanisms, more NS patients and their families should be investigated by means of NGS approach. This might be instrumental for the discovery of subclinical variants in functionally related genes that are apparently benign if inherited alone, but potential pathogenic if co-segregating in a single patient.

In conclusion, our results together with the growing body of literature, indicate that the routinely application of massive sequencing techniques is uncovering the emerging contribution of hypomorphic variants to the clinical phenotype, unlike the application of Sanger method that led to interrupt the analysis when an affecting function variant was detected. Geneticists need to be aware of the complexity of the genetic characterization and counseling about NS and RASopathies. Therefore, the revaluation of many monogenic Mendelian conditions characterized by a broad spectrum of phenotypes, such as NS, may lead to the identification of different inheritance models, able to explain variable expressivity and incomplete penetrance.