Novel Genetic and Biochemical Insights into the Spectrum of NEFL-Associated Phenotypes

Patients

Patients included in this study were phenotyped in the Department of Pediatric Neurology of the University Hospital of Essen (Duisburg-Essen University). Written informed consent for clinical description, genetic studies, and utilization of the muscle biopsy for research purposes was obtained from the patients’ parents. The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Ethics Committee of University of Duisburg-Essen (19-9011-BO).

Molecular genetic studies

Genomic DNA was obtained from total peripheral blood samples using routine extraction methods. Next-generation sequencing (gene panel analysis including 1.678 genes with relevance in clinical diagnostics, custom-made target enrichment, Agilent Sure Select) was carried out on an Illumina NextSeq 500 system (Illumina, San Diego, CA) as 150 bp paired-end runs using v2.0 SBS chemistry. Sequencing reads were aligned to the human reference genome (GRCh37/hg19) using BWA (v0.7. 13-r1126) with standard parameters. Statistics on coverage and sequencing depth on the clinical targeted regions (i.e. RefSeq coding exons and +/–5 intronic region) were calculated with a custom script. SNV and INDEL calling on the nuclear genes was conducted using SAMtools (v1.3.1) with subsequent coverage and quality dependent filter steps. Variant annotation was performed with snpEff (v4.2) and Alamut-Batch (v1.4.4). Only variants (SNVs/small INDELs) in the coding and flanking intronic regions (±50 bp) were evaluated. BayesDel (https://fenglab.chpc.utah.edu/BayesDel/BayesDel.html) and the thresholds published by Pejaver et al. [27] were used for bioinformatic predictions. Variants were classified according to the ACMG (American College of Medical Genetics and genomics) guidelines taking the current recommendations from ClinGen's Sequence Variant Interpretation Working Group (SVI) into account [16, 28] and classifying variants with the 5-tier classification system: class 5 (pathogenic), class 4 (likely pathogenic), class 3 (uncertain variants or variants of unknown significance, VUS), class 2 (likely not pathogenic) and class 1 (not pathogenic).

In silico modelling of amino acid substitutions affecting Pro 8

ColabFold version 1.5.2 was employed using the AlphaFold2_mmseqs2 Jupyter notebook to predict the three-dimensional structures of the wild-type NEFL protein (UniProt accession: P07196) and five missense NEFL protein variants [29]. Each three-dimensional structure was predicted alongside pLDDT score, which is a per-residue confidence metric between 1–100 generated by AlphaFold2 [30]. PyMOL was then used to individually superimpose the wild-type NEFL protein model onto each of the missense protein variant models using the “super” command. The “super” command also computes a root-mean-square deviation (rmsd) score for each superimposition, representing structural similarity between the two models [31]. Finally, the predicted protein models were uploaded to PremPS, a web server that can be used to predict the effects of missense mutations on protein stability, measured by unfolding free energy change (ΔΔG) [32]. The five NEFL missense mutations were specified to the PremPS server, and the changes in non-covalent interactions between wild-type NEFL and each missense variant were modelled, exclusively focusing on non-covalent interactions between the mutated amino acid residue and its adjacent residues. A ΔΔG value (kcal/mol) was estimated for each mutation as an indicator of protein stability change as a result of the mutation, with positive and negative ΔΔG values indicating mutations that either destabilize or stabilize the protein structure, respectively [32].

Biopsy work-up

Muscle biopsy specimen investigated in this study (from patient 1) has initially been collected for diagnostic purposes including histology, enzyme histochemistry, immunofluorescence and immunohistochemical investigations. Serial cryosections (7μm) of transversely-oriented muscle blocks were stained according to standard procedures with hematoxylin and eosin (H&E), Gömöri trichrome (GT), oil red, COX-SDH and SDH and nicotinamide adenine dinucleotide tetrazolium reductase (NADH-TR) [33]. Microscopy was performed using a Zeiss Axioplan epifluorescence microscope and a Zeiss Axio Cam ICc 1.

Glutaraldehyde-fixed specimens were processed for ultrastructural examination by standard procedures. The tissue was post-fixed in 1% osmium tetroxide and embedded in Epon 812. Semithin sections for light microscopy were stained with toluidine blue. Ultrathin sections were contrasted with uranyl acetate and lead citrate and examined using a Philips CM10 transmission electron microscope as described [34].

Coherent anti-Stokes Raman scattering microscopy (CARS)

From a cyro-preserved muscle biopsy specimen, 5μm sections were prepared and used for coherent anti-Stokes Raman scattering (CARS) as well as for second harmonic generation (SHG) microscopy. Samples were dried in a nitrogen stream and required no further sample preparation.

For the CARS/SHG measurements, a modified Leica TCS SP 8 CARS was used with an APE picoEmerald as the laser source (as described previously, [35]). CARS and SHG signals were detected simultaneously in backward (Epi; E) and in forward (F) direction. For the imaging, we used a 40× water immersion objective lens (IRAPO 40×/1.10 WATER). For imaging and acquisition of CARS spectra, a field of 291×291μm (2048×2048 pixels) was measured on the samples. By tuning the pump laser from 804.0 nm to 826.4 nm with a step size of 0.7 nm, a pixel dwell time of approximately 10μs and averaging of two images, we acquired the FCARS spectra.

Proteomic profiling

Proteomic profiling of whole muscle protein extracts was carried out as described previously in a data-independent-acquisition mode [36].

Clinical findings

Patient 1: was born as the first child of non-consanguineous parents of German origin. Pregnancy and birth were uneventful. At the age of 1.5 years, developmental language problems occurred, sensorineural hearing loss was diagnosed and linked to the compound-heterozygous CX26 (GJB2) variants c.109G > A, p.(Val37Ile) (ClinGen Hearing Loss Expert Panel Curation classified as “pathogenic”) and c.100_101delinsCC, p.(Met34Pro) (classified as “likely pathogenic” using ACMG-AMP codes PM5, PM3, PM1, PM2_SUP). However, motor milestones were achieved in time. Due to cognitive impairment (IQ: 64) and behavioral problems (both are not part of the phenotypic spectrum associated with variants in CX26) a basic diagnostic work-up (echocardiography, electrocardiogram [ECG], laboratory values including basic laboratory examinations for metabolic disorders as well as a screening for congenital disorders of glycosylation (CDG) was performed with normal results. Karyotype was 46, XY. Electroencephalography (EEG) at 24 months of age showed no pathological findings, brain and spine magnetic resonance imaging (MRI) revealed normal findings. At the age of 2 years, muscle cramps and muscular stiffness of lower extremities (especially in the morning) as the first signs of a neuromuscular involvement occurred. Creatine kinase (CK) was slightly elevated (210 U/l; norm < 149 U/l). Muscle biopsy was performed at the age of 6 years assuming a muscular disease. Phenotyping at his last visit (at the age of 11 years) revealed reduced endurance (walking distance < 300 m), muscle cramps under muscular strain accompanied by slowly progressing gait disturbances, hyperlordosis, normal reflex levels for patellar tendon reflex (PTR) and absent achilles tendon reflex (ATR), pes cavus, small CMAP, marginally reduced nerve conduction velocities (NCV) and prolonged distal latencies (Fig. 1).

Fig. 1

Clinical findings in one NEFL-patient. Patient 1 at the age of 3 years (A and D) and 11 years (B, C and E). Note lack of progress in pes cavus within the period of 7 years. Hyperlordosis at the age of 11 years (C).

Patient 2: was born as the first child of non-consanguineous parents of German origin. She was preterm born at 28 weeks of gestation due to cardiac failure of the mother under resuscitation. Mother deceased during delivery. After 10 weeks of hospitalization the infant was discharged with need of tube feeding for 3 months. She had intraventricular bleeding one sided due to prematurity (ultrasound). MRI of brain at the age of 11 years was normal. Persistent ductus arteriosus and atrial septum defect were operated during infancy. Her motor milestones were slightly delayed; she was able to walk unsupported at the age of 21 months, but falls were frequent. Speech and cognitive development were age appropriate. At the age of 5 years, the patient developed an abnormal gait pattern with unstable gait, with inward positions of the knee and the development of pes cavus. She had reduced proximal muscle strength and could stand up from the sitting position only with support. She had problems to open a bottle due to distal muscle weakness. This new symptom could not be explained as a consequence of premature birth and thus further diagnostic was initiated. The basic work-up (echocardiography, electrocardiogram [ECG], laboratory values including CK and basic laboratory examinations for metabolic disorders as well as electroencephalography [EEG], brain and spine MRI) revealed normal findings. At the age of 11 years, the patient underwent surgery for feet deformity including heel cord lengthening. Her clinical findings at her last visit (at the age of 16 years) revealed reduced walking distance (approximately 1200 m), brisk reflex levels, sensory ataxia, pes cavus, small CMAP, reduced NCV and prolonged distal latencies in her arms (Nervus (N) medianus and N. ulnaris, sensible and motoric) and absent in legs (motoric and sensible NCV not reproducible) (patient 7 in Table 1).

Father of patient 2: He reported pes cavus at both feet since early infancy/childhood (exact age of onset not known) without any progression, no reduced endurance or gait disturbances. His father had also pes cavus as far as he can remember. Last neurological examination at the age of 47 years showed normal results for his muscular strength and muscular tonus with brief muscle reflexes in his upper extremities and PTR, but reduced reflexes at ATR. Nerve conduction velocity studies revealed prolonged distal latencies (6.0/6.3 m/sec) and slightly reduced motor NCV (40/37 m/sec) at both N. peroneus as a mild neuropathic sign (patient 8 in Table 1).

Molecular genetic findings

Patient 1: By next generation sequencing we identified the heterozygous variant c.310T>G; p.(Phe104Val) in the NEFL gene (NM 006158.4). The nucleotide substitution from T to G at position 310 leads to a predicted amino acid exchange from phenylalanine to valine. The corresponding amino acids show minor differences in polarity, charge and size (Grantham dist.: 50 [0–215]). The amino acid phenylalanine at position 104 is located in the “intermediate filament protein” domain of the NEFL protein and is highly conserved across many species. The absence of this variant in population-specific databases can be considered as supporting evidence for a pathogenic effect. This variant has been detected in one individual from a cohort of patients with peripheral neuropathy [37]. Bioinformatic prediction using BayesDel suggest a strong pathogenic effect. The result of the family analysis suggests a de novo origin of the sequence variant in the index patient, supporting a pathogenic effect of the variant. According to the ACMG-guidelines [28], the variant was classified as “likely pathogenic” (evidence codes: PS4_SUP,, PM2_SUP, PM6_SUP, PP3_Strong).

Patient 2: We identified the heterozygous variant c.22C > T; p.(Pro8Ser) in the NEFL gene. This variant has not yet been described in the literature, but three other amino acid exchanges at the same codon have been described as being pathogenic (c.23C > G; p.(Pro8Arg); c.23C > A; p.(Pro8Gln) & c.23C > T; p.(Pro8Leu)) [38], indicating that the proline at this amino acid position is critical for NEFL protein activity. Assumed pathogenicity of the variant is not only suggested by the biochemical nature (substitution of a neutral-polar but non-circular amino acid (Ser) to a polar circular one (Pro)) but was also supported by segregation analysis: A targeted Sanger sequencing revealed that the variant c.22C > T; p.(Pro8Ser) could be detected in the mildly affected father as a low-grade somatic mosaic (about 15% in peripheral blood). Additionally, this variant has been detected in a confirmed de novo constellation in a Taiwanese CMT-patient [39] and has moreover been described in at least two unrelated individuals with NEFL-associated phenotypes (ClinVar VCV000576990.5). However, in case of this one reported individual, the authors do not provide any data supporting pathogenicity and therefore we have at least included this “evidence” under PS4 (proband counting). Taken together, this variant was classified as “likely pathogenic” (evidence codes: PS2_MOD, PS4_SUP,, PM2_SUP, PM5), according to the ACMG guidelines and following the recent recommendation from ClinGen's Sequence Variante Interpretation Working Group (https://clinicalgenome.org/working-groups/sequence-variant-interpretation/).

In silico modelling of amino acid substitutions affecting proline at position 8

To address the impact of the nature of amino acid substitutions on the manifestation of clinical symptoms in terms of severity of phenotypes associated with different missense variants affecting amino acid position 8, Alphafold2-based in silico studies were carried out: the six NEFL protein structures were predicted with pLDDT scores ranging from 71.8–72.4 (Fig. 5A). When superimposed onto the wild-type NEFL protein model, the Phe104Val and Pro8Leu protein variants exhibited rmsd scores of 4.739 Å and 7.959 Å, respectively, which were the highest rmsd scores among the modelled variants (Fig. 5B & Table 1). On the other hand, Pro8Ser, Pro8Gln, and Pro8Arg exhibited rmsd scores between 2.5 Å and 3.0 Å when superimposed onto wild-type NEFL (Fig. 5B & Table 1). Identical structures are indicated by a rmsd of 0 Å, a rmsd up to 2.0 Å reflects good structural similarity, and a rmsd of 3.0 Å and above indicates low similarity [40–42]

PremPS predicted that all five variants alter the hydrophobic interactions between the mutated and its adjacent residues (Fig. 5C). The p.(Phe104Val) variant introduces a new hydrophobic interaction between Val104 and Val108. The Pro8Ser variant lacks all three hydrophobic interactions between Pro8 and Tyr6 that are present in wild-type NEFL. Pro8Arg leads to the loss of a hydrophobic interaction between Arg8 and Tyr6, and p.(Pro8Gln) leads to the loss of a hydrophobic interaction between Gln8 and Tyr6. Lastly, the p.(Pro8Leu) mutation causes one hydrophobic interaction to be lost between Leu8 and Tyr6 while also introducing two new hydrophobic interactions between Leu8 and Tyr10. All changes in hydrophobic interactions were associated with a ΔΔG value (kcal/mol) (Table 1). PremPS estimated either stabilising or destabilising ΔΔG values for the five NEFL variants. Both stabilising and destabilising ΔΔG values have been associated with deleterious mutations. The Pro8Leu variant showed a ΔΔG of –0.35 kcal/mol compared to wild-type NEFL, which was the highest magnitude ΔΔG value among the variants tested. Thus, our in silico predictions underline the pathogenic character of the amino acid substitution identified in our NEFL-patients and moreover demonstrate that missense variants affecting the same amino acid substitution may have diverging deleterious effects in terms of irregular stabilization or destabilization of NEFL.

Histological and electron microscopic findings on NEFL-mutant muscle

Histological examination did not reveal striking myopathic changes in the patient 1’s muscle biopsy (H&E staining shown in Fig. 2A) except the indication of altered mitochondrial distribution based on NADH-TR staining (Fig. 2B). However, ultra-structural studies utilizing electron microscopy unraveled focal alterations of sarcomere pattern with streaming of Z-discs and incidental pronounced Z-disc disruption as minicore-like lesions in addition to small cytoplasmic bodies (Fig. 2C & D). No larger autophagic vacuoles, or tubular aggregates were identified. Mitochondria presented with normal size and shape and no ultrastructural abnormalities of nuclei were found (data not shown).

Fig. 2

Microscopic findings in NEFL-mutant muscle. Haematoxylin & Eosin (H&E) staining did not reveal striking myopathic changes (A) whereas NADH-TR in some muscle fibers revealed increased subsarcolemmal reactivity indicative for mitochondrial accumulations in addition to pale regions in sum indicating altered mitochondrial distribution (B, arrows). Electron microscopic studies revealed focal alterations of sarcomere pattern with streaming of Z-discs and incidental pronounced Z-disc disruption as minicore-like lesions and small cytoplasmic bodies (C and D, arrowhead).

Coherent anti-Stokes Raman scattering microscopy (CARS)

CARS offers the unique opportunity to investigate protein and lipid imbalances in an unbiased manner by making use of a minimal amount of material and is thus a valuable analytical technique in the investigation of pathological features in biopsy material such as muscle biopsy specimens. Thus, we applied this label-free technique to obtain further unbiased insights into the muscle cell vulnerability underlying the pathogenic NEFL-variant. To locate spectroscopic features in the muscle biopsy samples (NEFL-patient 1 and controls), first the spectrum of the fibre background (areas in which no features) was defined, and a mean spectrum was calculated. The fibre background of NEFL-patient 1 (Fig. 3A, pBG) shows a lower signal intensity at 2921 cm^–1 than at 2868 cm^–1. The signal at 2921 cm^–1 is characteristic for proteinogenic content [43, 44]. The standard peak of ordered lipids [45] at 2889 cm^–1 seems to be shifted to 2868 cm^–1. In contrast, in subsarcolemmal regions of NEFL-mutant muscle fibres, the intensity at 2921 cm^–1 and 2868 cm^–1 is comparable (Fig. 3A, pProt). For CARS A, a higher intensity for 2921 cm^–1 based spectra can be observed in the subsarcolemmal regions compared to the fibre background. In the light of the above-mentioned aspects, we concluded based on these data that in NEFL-mutant muscle demonstrates an increased protein content (mirrored by increased signal at 2921 cm^–1) is present in some regions of fibers.

Fig. 3

Spectroscopic findings in NEFL-mutant muscle obtained by F-CARS measurements. Protein-containing regions in muscle fibres of NEFL-patient 1 at 2921 cm^–1 (marked with arrows; scale bar 10μm) (A). Subsarcolemmal regions in muscle fibres of disease control present an increased signal at 2921 cm^–1 (marked with arrows; scale bar 15μm) (B). Lipid features identified within and between muscle fibres of disease control with increased signal at 2847 cm^–1 (marked with arrows; scale bar 15μm). Within the fibres, the lipid is present as sarcoplasmic dots (C). Muscle fibres of a representative healthy control (scale bar 10μm). Areas representing the fibre background without features (pBG) and the protein-containing subsarcolemmal regions (pProt) were selected from the measurements carried on the NEFL patient muscle biopsy specimen. Spectra were collected from these areas. Likewise, spectra were taken at disease control from the fibre background (dcBG), protein-containing subsarcolemmal region (dcProt), and lipid-containing region (dcLip). For each group, spectra were normalized to 1 and averaged. The standard deviation is shown in grey. The vertical lines mark the wavenumbers 2847 cm^–1, 2868 cm^–1, and 2921 cm^–1 (D). Comparison of the average muscle fibre calibres. NEFL-patient 1 : 24.01μm±4.94μm (721 fibres); disease control: 16.97μm±9.78μm (782 fibers); controls: 1 : 39.38μm±7.23μm (318 fibers); 2 : 43.48μm±10.49μm (184 fibers); 3 : 72.82μm±16.29μm (436 fibers); 4 : 19.62μm±3.62μm (430 fibers); 5 : 36.94μm±6.09μm (89 fibers); all controls: 43.92μm±23.29μm (1457 fibers) (E).

Spectral analysis of the disease control (neurogenic muscular atrophy in the biopsy of an age-matched child) revealed differences in comparison to the NEFL-patient 1 (Fig. 3B & C). The fibre background (dcBG) shows a slightly higher intensity at 2921 cm^–1 than at 2868 cm^–1. Thus, the spectrum is similar to that of the controls (Fig. 3D, cBG). Also, less pronounced subsarcolemmal spots with significantly increased signal at 2921 cm^–1 were found in the disease control (Fig. 3, dcProt), but regions in and between fibres were also identified with high signal at 2847 cm^–1 (Fig. 3C, dcLip). Especially disordered lipids show a strong signal at 2847 cm^–1. No lipid features were found, neither in controls nor in the NEFL-patient.

To investigate muscle fibre calibre, we measured the average of muscle cells in the biopsy derived from NEFL-patient 1 in comparison to a disease control (denervation atrophy) and normal controls: while the averaged muscle fibre calibre in the NEFL-patient 1’s muscle biopsy is 24.01μm±4.94μm (721 fibres), a calibre of 54.09μm±20.33μm (1027 fibres) was detected in the normal controls and 16.97μm±9.78μm (782 fibres) in the disease control (Fig. 3E). Thus, our CARS-based investigations also enabled to demonstrate smaller muscle fibre calibre in NEFL-patient 1, a microscopic finding which accords with the assumption of neuropathy-related denervation.

Proteomic profiling of NEFL-mutant muscle derived from patient 1

To identify the proteins (reflecting molecular processes) underlying the pathogenic NEFL variant in skeletal muscle, we performed label-free untargeted proteomic profiling on quadriceps muscles derived from patient 1 in comparison to three healthy age- and gender-matched controls (Fig. 4A). By measuring 5777 unique peptides (Fig. 4B), we identified 58 proteins with an increased abundance of more than 1.3-fold whereas 60 proteins were identified with a decrease of a minimum of 0.7-fold (Fig. 4C & Table 2). These 118 altered proteins were quantified based on at least two unique peptides and significantly dysregulated with a p-ANOVA≤0.05 (Table 2). All functional information about these dysregulated proteins listed in Table 1 was extracted from uniprot (www.uniprot.org, accessed on 10th of March 2023). We performed GO term analysis focussing on the up- and down-regulated proteins separately (Fig. 4D): within the group of upregulated proteins, the most altered biological processes include organization of the extracellular matrix including collagen fibril organization (indicative for fibrotic processes), activity of complex I of mitochondrial respiratory chain and metabolic processes of reactive oxygen species and carbohydrates, tricarboxylic acid cycle, axogenesis including axon extension involved in axon guidance, chaperone-mediated protein folding and muscle contraction (Fig. 4B). Molecular functions affected by increased proteins contain oxidoreductase activity including NADH dehydrogenase activity and electron carrier activity, extracellular matrix structural constituent and matrix binding, double-stranded RNA binding, activities of superoxide dismutase and thioester hydrolase, binding of proteins to chaperones (to warrant proper folding) and muscle alpha-actinin protein binding (Fig. 4D). In accordance with affected biological processes and molecular functions, affected cellular compartments include the extracellular matrix, mitochondria, lumen of Golgi apparatus and lysosomes, peroxisomes, nuclear chromatin, sarcolemma including membrane rafts, actin cytoskeleton and the myelin sheath (Fig. 4D). A GO-term-based pathway analysis of down regulated proteins revealed that a significant proportion of these proteins affect different biological processes including glycolysis, tricarboxylic acid cycle, mitochondrial acetyl-CoA biosynthetic processes, creatine metabolic processes, muscle organ development and functions such as muscle contraction and myofibril assembly. Molecular functions affected by decreased protein abundances contain structural constituent of muscle and cytoskeleton (including binding of cytoskeletal proteins such as titin and tropomyosin), ion channel binding as well as activities of structural proteins and pyruvate dehydrogenases (Fig. 4D). In line with the affected biological processes and molecular functions decreased proteins impact on myelin sheath, extracellular matrix, myonuclei, mitochondria and the sarcomeric system (Z-disc, I and M bands, actin and thin filaments) (Fig. 4D).

Fig. 4

Proteomic profiling of NEFL-mutant muscle. (A) Schematic representation of our workflow applied to decipher the proteinogenic changes in NEFL-mutant muscle. (B) Statistics of our proteomic findings. (C) Volcano-plot indicating the statistically significant dysregulated proteins. Increased proteins are highlighted as purple dots whereas decreased proteins are depicted as orange dots. (D) Results of in silico based (GO-term) analysis of the proteomic signature of NEFL-mutant muscle indicating biological processes, molecular functions and subcellular structures affected by the dysregulated proteins, respectively.

Fig. 5

In silico studies on amino acid substitutions affecting NEFL. (A) ColabFold predictions of three-dimensional protein structure of (i) Wild-type (WT) NEFL (ii) Phe104Val (iii) Pro8Ser (iv) Pro8Arg (v) Pro8Gln (vi) Pro8Leu. Residues are colour-coded on a spectrum of dark blue to red based on pLDDT score, with blue representing high pLDDT and red representing low pLDDT. WT NEFL pLDDT: 72.2, Phe104Val pLDDT: 72.2, Pro8Ser pLDDT: 72.4, Pro8Arg pLDDT: 72.2, Pro8Gln pLDDT: 71.8, Pro8Leu pLDDT: 72. (B) Root-mean-square deviation (Rmsd) scores calculated from superimpositions of each NEFL protein variant with the WT NEFL protein structure. Superimpositions were performed in PyMOL. (C) PremPS predictions of differences in hydrophobic interactions between mutated residues and adjacent residues in NEFL variants compared to WT NEFL.