Investigating genotype–phenotype relationship of extreme neuropathic pain disorders in a UK national cohort Open Access

Abstract

The aims of our study were to use whole genome sequencing in a cross-sectional cohort of patients to identify new variants in genes implicated in neuropathic pain, to determine the prevalence of known pathogenic variants and to understand the relationship between pathogenic variants and clinical presentation. Patients with extreme neuropathic pain phenotypes (both sensory loss and gain) were recruited from secondary care clinics in the UK and underwent whole genome sequencing as part of the National Institute for Health and Care Research Bioresource Rare Diseases project. A multidisciplinary team assessed the pathogenicity of rare variants in genes previously known to cause neuropathic pain disorders and exploratory analysis of research candidate genes was completed. Association testing for genes carrying rare variants was completed using the gene-wise approach of the combined burden and variance-component test SKAT-O. Patch clamp analysis was performed on transfected HEK293T cells for research candidate variants of genes encoding ion channels. The results include the following: (i) Medically actionable variants were found in 12% of study participants (205 recruited), including known pathogenic variants: SCN9A(ENST00000409672.1): c.2544T>C, p.Ile848Thr that causes inherited erythromelalgia, and SPTLC1(ENST00000262554.2):c.340T>G, p.Cys133Tr variant that causes hereditary sensory neuropathy type-1. (ii) Clinically relevant variants were most common in voltage-gated sodium channels (Na_v). (iii) SCN9A(ENST00000409672.1):c.554G>A, pArg185His variant was more common in non-freezing cold injury participants than controls and causes a gain of function of Na_V1.7 after cooling (the environmental trigger for non-freezing cold injury). (iv) Rare variant association testing showed a significant difference in distribution for genes NGF, KIF1A, SCN8A, TRPM8, KIF1A, TRPA1 and the regulatory regions of genes SCN11A, FLVCR1, KIF1A and SCN9A between European participants with neuropathic pain and controls. (v) The TRPA1(ENST00000262209.4):c.515C>T, p.Ala172Val variant identified in participants with episodic somatic pain disorder demonstrated gain-of-channel function to agonist stimulation. Whole genome sequencing identified clinically relevant variants in over 10% of participants with extreme neuropathic pain phenotypes. The majority of these variants were found in ion channels. Combining genetic analysis with functional validation can lead to a better understanding as to how rare variants in ion channels lead to sensory neuron hyper-excitability, and how cold, as an environmental trigger, interacts with the gain-of-function Na_V1.7 p.Arg185His variant. Our findings highlight the role of ion channel variants in the pathogenesis of extreme neuropathic pain disorders, likely mediated through changes in sensory neuron excitability and interaction with environmental triggers.

Graphical abstract

Open in new tabDownload slide

Introduction

Neuropathic pain occurs as a consequence of a disease or lesion in the somatosensory nervous system.¹ It affects 6.9–10% of the general population² and has a harmful impact on physical health, psychological health and quality of life.³ Understanding the role of genetic factors in neuropathic pain may reveal new pathophysiological mechanisms and is under-explored.^4,5

Extreme pain phenotypes, caused by rare high-impact genetic mutations, offer insight into fundamental neurobiological mechanisms of pain.⁶ The phenotypes can range from congenital insensitivity to pain⁷ to enhanced pain perception. Different types of mutations in the same gene can cause a spectrum of phenotypes. For example, bi-allelic loss of function mutations in SCN9A, the gene encoding the voltage-gated sodium channel (Na_v) 1.7 and which is highly expressed in peripheral sensory neurons,⁸ causes congenital insensitivity to pain.⁹ In contrast, monoallelic gain-of-function variants in the same gene is associated with pain disorders, which are inherited in a Mendelian fashion. These include inherited erythromelalgia¹⁰ and paroxysmal extreme pain disorder.¹¹

Variants in genes causing Mendelian pain disorders may act as risk factors for common acquired neuropathic pain disorders. For example, SCN9A variants are implicated in more common neuropathic pain disorders such as idiopathic small fibre neuropathy¹² and painful diabetic neuropathy.¹³ The SCN9A NM_002977.3:c.3448C>T, p.Arg1150Trp variant modulates risk and severity of pain across different chronic pain disorders.¹⁴ An environmental trigger that interacts with genes may cause neuropathic pain as some variants are common in the general population. Identification of such gene variants is important in diagnosis, genetic counselling and treatment enabling a stratified approach, such as the use of sodium channel blockers (e.g. lacosamide) in patients with small fibre neuropathy and gain-of-function Na_v 1.7 variants.¹⁵ In some cases, such as inherited erythromelalgia, a personalized management approach can be used.^16,17 Furthermore, there are rare inherited conditions where direct treatment can arrest progression, such as Fabry’s disease and hereditary transthyretin amyloidosis. In these diseases neuropathic pain is often the first symptom, the disease will progress if untreated and timely genetic diagnosis is essential to initiate appropriate treatment.¹⁸

The National Institute for Health and Care Research (NIHR) BioResource Rare Disease project applied whole genome sequencing (WGS) to a range of rare diseases, including neuropathic pain disorders.¹⁹ The aims of our study were to aid the genetic diagnosis of patients with extreme neuropathic pain phenotypes, to determine the prevalence of gene variants associated with neuropathic pain disorders and to understand how the functional changes caused by gene variants relate to clinical presentation and neuropathic pain (Supplementary Fig. 1).

Methods

Recruitment and clinical phenotyping of participants

We recruited patients with extreme neuropathic pain phenotypes, both sensory loss and gain, from secondary care clinics in the UK, located in Oxford, London, Salford and Newcastle. Study participants with a history of lifestyle-altering sensory disorder, either pain or loss of sensation, for greater than 3 months were invited to participate. The criteria for clinical case definitions are shown in Table 1.¹⁹ We excluded patients with a known underlying genetic cause of chronic pain, e.g. Fabry’s disease and SCN9A congenital erythromelalgia (genetic pre-screening for these disorders was not mandatory), pregnancy, coincident major psychiatric disorders, poor or no English language skills, patients with documented central nervous system lesions, or patients with insufficient mental capacity to provide informed consent or to complete phenotyping. Description of the clinical phenotyping can be found in the summary NIHR Bioresource paper by Turro et al.¹⁹ and is briefly described here.

Study participants attended a single appointment that included a clinical assessment, screening for neuropathic pain and specialized investigations to investigate distal symmetrical polyneuropathy. A detailed medical and drug history was taken, followed by a structured upper and lower limb neurological examination to detect clinical signs of distal symmetrical polyneuropathy.^20,21 DN4 questionnaire was used as a screening tool for neuropathic pain.²² Confirmatory tests included nerve conduction studies,²³ skin biopsy for intra-epidermal nerve fibre density^24,25 and thermal thresholds in the area of neuropathic pain.²⁶ Study participants’ pain was assessed and graded (Supplementary Fig. 2) according to published guidelines.²⁷

Whole-blood samples were collected and sent to the NIHR BioResource laboratory in Cambridge. A detailed description of the DNA sequencing, WGS data-processing pipeline and identification of relevant gene variants can be found in Turro et al.¹⁹ and relevant aspects are summarized below.

All participants provided written informed consent in accordance with the Declaration of Helsinki. The study was approved by the East of England Cambridge South national research ethics committee (REC) reference 13/EE/0325.

Analysis plan

Genetic analysis was completed in two parts. The first analysis was to identify variants of clinical relevance in known pain genes., grouping rare variants in genes for all neuropathic pain phenotypes considering both a targeted panel of pain genes and their promoters and all genes in the human genome that carried rare variants. Two ion channel variants were selected for electrophysiological analysis to investigate their functional impact (Supplementary Fig. 1).

Clinical reporting of pertinent findings

Gene list and transcript selection

A list of genes separated into three tiers were curated at the time of recruitment in 2015 (Tables 2–4) and updated in 2021 (Table 5). The division was based on the strength of evidence for the gene being linked to neuropathic pain.¹⁹ Only Tier 1 genes (Table 2) were discussed in the multidisciplinary team meetings and considered for clinical reporting.

Variant filtering to identify variants of clinical relevance

Variants of the 14 Tier 1 neuropathic pain disorders genes were prioritized based on

• Minor Allele Frequency (MAF) in gnomAD < 1/1000 if the variant was not previously described in association with disease, or<25/1000 if the variant was present as disease-causing or questionable disease-causing in Human Gene Mutation Database;

• predicted impact according to Ensembl Variant Effect Predictor were ‘High’, ‘Moderate’ or ‘splice region variant’. Variants with more than three alternate alleles or an internal MAF > 10% were discarded to guard against errors in repetitive regions and prevent potential systematic artefacts.¹⁹

Variant interpretation in multidisciplinary teams meeting

Prioritized variants were assessed by a multidisciplinary team.¹⁹ Pathogenicity assignment was ascribed according to guidelines of The American College of Medical Genetics.²⁸ Variants were classified as pathogenic, likely pathogenic, variant of uncertain significance (VUS), likely benign or benign. Clinically relevant variants were those deemed medically actionable, i.e. variants that could result in specific, defined medical recommendations, and were reported to the referring clinician.

Pathogenic and likely pathogenic variants were deemed clinically relevant.

VUS were deemed clinically relevant if

1. in vitro functional studies showed that the variant altered function but the relationship to clinical phenotype was not clear, particularly if a significant gene and environment interaction was suspected; or

2. in silico analysis of the variant suggested a significant effect on function as defined by

   • variant position: in a functionally important domain; located in a highly evolutionary conserved section; significant biochemical consequences of amino acid exchange; affecting protein structure;

   • minor allele frequency;

   • concordant pathogenicity through multiple in silico algorithms such as SIFT, Polyphen and Align GDVD.

VUS that did not meet the above criteria, likely benign and benign variants were not reported.

Group-wise rare variant association testing

Genetic association testing in genes carrying rare variants was carried out using the gene-wise approach of the combined burden and variance-component test SKAT-O.^29,30 The SKAT-O test combines a standard gene burden test that maximizes power under the assumption that all rare variants collapsed in a specific gene region are causal and acting in the same direction towards the phenotype and the sequence kernel association that calculates the weighted sum of squares of the variant score statistics and thus is more robust to the presence of variants with opposing effects. The combined SKAT-O test calculates a linear combination of the burden and variant-component test. Parametric bootstrap was used to resample residuals under the null model. The same resampled bootstrap phenotypes, considering covariates, were used for each gene. The weighting of the linear combination is optimized from the data itself. The Rho statistic indicates the weighting of each test and gives an indication of the causality and directionality of the variants’ effects, with rho = 1 reducing SKAT-O to a burden test (high percentage of causality in the same direction) and a rho = 0 to a SKAT (causal and non-causal variants with opposing directions). Age, sex and the three first Principal Components of genetic variation were used as covariates. Only unrelated individuals were used in the group-wise rare variant analysis. Analysis was done per ethnicity and neuropathic pain clinical phenotype.

Data filtering and pre-process

Only participants with neuropathic pain were included for gene-wise rare variant association testing. In total 39 participants were excluded from the original cohort of 205 for the SKAT-O analysis (8 with no neuropathic pain; 2 age not available; 7 sequencing by Genomics England not available as mapped to Grch38 genome assembly and joint genotyping not conducted; 4 as only founders from families included). Exclusions were also based on ethnicity. We only considered 128 Europeans and 38 Africans with neuropathic pain as cases, 4 South Asians and 14 people of ‘Other’ ethnicity were not included. Unrelated individuals from the NIHR Bioresource Rare Diseases cohort, not recruited as a part of neuropathic pain disorders nor the Neurodevelopmental Disorders cohort, were included as controls. Analysis for individuals with European and African ethnic origin was separated. The sample size used was 128 Europeans and 38 Africans as cases versus 5945 and 154 as ethnically matched controls, respectively.

VCF files were normalized and left aligned, and multi-allelic single nucleotide polymorphisms (SNPs) were broken down to bi-allelic using BCF-tools. Only SNPs and indels that were rare in gnomAD (MAF < 0.001), not common in the whole NIHR Bioresources Rare Diseases cohort (MAF < 0.05), and had MAF < 0.05 in cases and controls combined for each phenotype were considered. Variants should have passed the quality filters in the joint genotyping calls and be of high quality with a call rate > 99%, following HWE equilibrium (P > 0.05). Alleles should have a high average base call depth (>10) and average genotype quality (>20). Regulatory regions for Tier 1–3 genes were downloaded using the GeneHancer resource form GeneCards. We then selected regions annotated as promoters or promoters/enhancers for the Tier 1–3 genes and collapsed in one group associated with the respective gene. For the panel of Tier 1–3 genes and their promoters, we considered variants with ‘modifier’, ‘moderate’ and ‘high’ impact effects. We further selected variants that were Nonsynonymous or have the EPACTS functional annotations: Essential_Splice_Site, Normal_Splice_Site, Start_Loss, Stop_Loss, Stop_Gain, 5′ UTR and 3′ UTR. For the whole gene set (all genes in the human genome), we only considered variants of ‘moderate’ or ‘high’ Variant Effect Predictor impact on protein-coding genes. Manta and Canvas software packages were used to detect deletions of >50 bp.¹⁹ Due to the limited sample size and multiple sub-phenotypes, we focused on protein-coding genes and high- and moderate-impact variants for the whole gene set analysis.

Separate analysis was completed for 36 Tier 1–3 genes (72 groups with their respective promoters); functionally validated variants in 3 voltage-gated sodium channels SCN9A, SCN10A and SCN11A⁸; and the whole gene set (∼20 000 genes).

Functional in vitro studies of SCN9A p.Arg185His and TRPA1 p.Ala172Val variants

Plasmids and site-directed mutagenesis

Human Na_V1.7 cDNA was cloned into a modified pcDNA3 expression vector containing downstream IRES and dsRED2 sequences (SCN9A-IRES-DsRED). Human β1 and β2 subunits were cloned into pIRES2-AcGFP (SCN1B-IRES-SCN2B-IRES-eGFP).⁹ Human TRPA1 cDNA was cloned into a modified pcDNA3 expression vector containing downstream IRES and dsRED2 sequences (TRPA1-IRES-DsRED).⁹ The mutations p.Arg185His, p.Ala172Val and p.Asn855Ser were introduced using QuickChange II XL site-directed mutagenesis kit (Agilent). The clones were sequenced by standard methods.

HEK293T cell culture and transfection

Human embryonic kidney HEK-293 T cells were grown in a Dulbecco’s modified Eagle’s culture medium (DMEM/F-12, Invitrogen) containing 10% foetal bovine serum and maintained under standard conditions at 37°C in a humidified atmosphere containing 5% CO₂. For the study of the SCN9A(ENST00000409672.1):c.554G>A, p.Arg185Hisvariant, cells were transfected using the jetPEI™ transfection reagent (Polyplus-transfection Inc.) with either WT or mutant Na_v1.7 channel combined with β1 and β2 subunits (2:1 ratio). For the study of TRPA1(ENST00000262209.4):c.515C>T, p.Ala172Val and TRPA1(NM_007332.3):c.2564A>G, p.Asn855Ser variants, cells were co-transfected with pMaxGFP (Amaxa) and either human TRPA1 wild type (WT) or human TRPA1- p.Ala172Val/p.Asn855Ser at a ratio of 1:5 to facilitate visualization of positively transfected cells. The total amount of plasmid DNA transfected was 1.2 μg/μl per 35 mm dish. pMaxGFP positive cells were used as control. Cells were used 36–72 h after transfection. Experiments were performed at room temperature and repeated on three or more separate transfections.

Electrophysiology

Voltage clamp experiments were performed on transfected HEK293T cells. Whole-cell patch clamp recordings were conducted using an Axopatch 200B Amplifier, the Digidata 1550B Low Noise Data Acquisition System and the pClamp10.6 software (Molecular Devices). Data were filtered at 5kHz and digitized at 20kHz. Capacity transients were cancelled and series resistance compensated at 70–90% in all experiments. Cells were continuously superfused with an extracellular solution or agonist-containing solutions through a common outlet. Standard extracellular solutions for the patch clamp experiments were used to study the respective SCN9A¹³ and TRPA1^31,32 variants.

Electrophysiology study of SCN9A p.Arg185His

The extracellular solutions contained (in mM): 140 NaCl,3 KCl, 1 CaCl₂, 1 MgCl₂, 10 HEPES, pH 7.3 with NaOH (adjusted to 320 mOsm/l with glucose). Patch pipettes were filled with an internal solution containing (in mM) 140 CsF, 10 NaCl, 1 EGTA, 10 HEPES, pH 7.3 with CsOH (adjusted to 310 mOsm/l with glucose) and had a typical resistance of 2–3 MΩ. Leak currents were subtracted using a P/5 protocol, applied after the test pulse. A holding potential of −100 mV and an intersweep interval of 10 s was used for all the protocols. Measurements were done at 10°C, 20°C and 30°C.

Electrophysiology studies of TRPA1 p.Ala172Val and p.Asn855Ser variants

To study the effects of agonist desensitization on channel activation, the extracellular solution contained (in mM): 127 NaCl, 3 KCl, 1 MgCl₂, 10 HEPES, 2.5 CaCl₂ and 10 glucose, pH 7.4 with NaOH. Osmolarity was adjusted to 310 mOsm/l with glucose. The intracellular solution contained (in mM): 135 KCl, 2 MgCl₂, 2 MgATP, 5 EGTA and 10 HEPES, pH 7.4 with CsOH. Osmolarity was adjusted to 300 mOsm/l with glucose.³¹

To study the effects of intracellular calcium in channel modulation, the extracellular solution contained (in mM): 140 NaCl, 4 KCl, 2 CaCl₂, 1 MgCl₂, 10 HEPES, pH 7.4 with NaOH. Osmolarity was adjusted to 310 mOsm/l with glucose. The intracellular solution contained (in mM): 130 KCl, 8 NaCl, 2 EGTA, 1 MgCl₂, 1 CaCl₂, 4 MgATP, 0.4 Na2GTP, pH 7.4 with CsOH. Osmolarity was adjusted 300 mOsm/l with glucose.

To investigate voltage dependence, currents were recorded during a voltage-step protocol consisting of 400 ms voltage steps to test potentials ranging from −100 to +180 mV, followed by a final invariant step to −75 mV (400 ms) to measure tail currents. The holding potential was set at −0 mV.

The voltage-dependence activation of hTRPA1 p.Ala172Val in response to mustard oil [Allyl isothiocyanate (AITC) a TRPA1 electrophilic agonist] and Menthol (non-electrophilic agonist of TRPA1) was measured. These recordings were performed in a calcium-containing extracellular solution, to preserve agonist desensitization properties.³³ Perfusion with TRPA1 agonists was performed through a custom-made gravity perfusion system. Immediately after establishing the whole-cell configuration, perfusion was switched to extracellular solution for 2 min before beginning the voltage clamp recording. AITC (Sigma 377430) was dissolved in DMSO (Sigma D2650), and Menthol (Sigma M2772) in ethanol. The maximum final concentration of either DMSO or ethanol did not exceed 0.1%. The effect of TRPA1 agonists on current–voltage curves was measured with a two-voltage-step protocol, as described above. Voltage ramps ranging from −100 to +100 mV for 500 ms, every 5 s, were applied to elucidate the temporal activation of hTRPA1 p.Ala172Val in response to AITC. In this case, the holding potential was set at −70 mV.

Current–voltage curves (I–V curves) were fitted using a combined Boltzmann and linear ohmic relationship: I/I_max = G_max (V_m−E_rev)/(1 + exp^{(V1/2−Vm)/k}). Normalized conductance–voltage curves (activation curves) were fitted with a Boltzmann equation G/G_max = 1/(1 + exp^{(V1/2−Vm)/K}), where G was calculated as follows G = I/(V_m−E_rev). Steady-state fast inactivation curves were fitted with I_T/I_Tmax = 1/(1 + exp^{−(V1/2−Vm)/k}). Tail current-derived voltage-activation curves were fitted to the Boltzmann equation: I_T/I_{T (Max)} = 1/(1 + exp^{[(Vm−V1⁄2)/k]}). In all the equations, V_1/2 represents the half-activation and half-inactivation membrane potentials; V_m is the membrane potential, E_rev the reversal potential, k the slope factor, G the conductance and I_T the current at a given V_m; G_max and I_Tmax are the maximum conductance and current, respectively; R_in is the fraction of channels that are resistant to slow inactivation. Statistical significance was set at P = 0.05 for group comparisons.

Statistical analysis

Electrophysiology data are presented as mean ± SEM. Statistical analysis for group comparisons included two-way ANOVA (temperature and genotype as a categorical variables) for p.Arg185His and one-way ANOVA (genotype as categorical variable) for p.Ala172Val with Sidak’s multiple comparison test. Statistically significant differences were defined as P < 0.05. Frequencies of individual genetic variants were compared across groups using a two-tailed Fisher’s exact test. Significance for gene-wise associations was set to the Bonferroni-adjusted 0.01 threshold for the number of genes considered. Unadjusted P-values are reported alongside the significance threshold.

Results

Study participants

A total of 205 study participants with extreme phenotypes were included (Fig. 1). Age ranged from 4.8 to 84.3 years, and 115 (56.1%) participants were men and 90 (43.9%) were women. In total, 38 participants were recruited of African descent. The majority of participants (90.2%) satisfied the criteria for probable or definite chronic neuropathic pain. Participants with possible neuropathic pain (5.9%) included those with ‘neuropathic-type’ pain (burning, stabbing, electric-like shocks, dysaesthesias) in a neuroanatomically plausible distribution, but with no evidence of nerve injury on clinical examination or specialized investigations. Examples include those with episodic pain syndromes or those with burning pain of the hands and feet with a normal clinical examination and investigations. A group of participants (8, 3.9%) did not experience neuropathic pain, including participants with loss of pain sensation.

Gene variants reported—Tier 1 gene analysis

After multidisciplinary team discussion, 26 (12.0%) gene variants were categorized as clinically relevant: 3 pathogenic (1.4%), 2 likely pathogenic (0.9%) and 21 VUS deemed relevant for reporting (9.7%). Apart from three participants who declined consent for feedback on genetic testing, all variants were reported to the referring clinician. The results are summarized in Table 6 and Supplementary Table 1.

Three participants were diagnosed with pathogenic mutations. Two sisters with severe erythromelalgia, present since childhood, have the same pathogenic variant in SCN9A p.Ile848Thr. This variant is not present in control populations, is reported in several inherited erythromelalgia pedigrees and causes gain of function through a hyperpolarizing shift in the voltage dependence of activation.^10,34 A pathogenic variant in SPTLC1 p.Cys133Trp was identified in a participant diagnosed with a painful sensorimotor neuropathy and is the commonest pathogenic SPTLC1 mutation causing hereditary sensory neuropathy type-1 identified in UK patients.³⁵

Two participants diagnosed with small fibre neuropathy carried likely pathogenic variants, SCN10A p.Gly1662Ser and SCN11A p.Cys1543Tyr. This conclusion was based on criteria available at the time of interpretation, low frequency in genetic databases, location in the important region of the channel, high probability of affecting channel function and previous description in small fibre neuropathy¹²; however, we note that this variant is classified as likely benign in ClinVar. The remaining 21 variants reported were classified as VUS, due to our conservative approach in assigning pathogenicity.

Of the 40 participants recruited with non-freezing cold injury, six participants (7.5% allele frequency) carry the SCN9A p.Arg185His variant, which is associated with small fibre neuropathy.^12,36 The six participants were all of the African descent (the majority from Ghana). Chronic non-freezing cold injury is an acquired painful sensory neuropathy observed almost exclusively in soldiers.³⁷ Soldiers of African descent are disproportionately affected when compared to Caucasian soldiers. The allele frequency of p.Arg185His variant is <1% in gnomAD (1.4% in the African/African-American population); 165 individuals carried the allele in the general gnomAD population, 125 of which were Africans (166 allele counts, 1 was homozygote). Within NIHR Bioresource controls of African ancestry, p.Arg185His allele frequency was 0.0065 (2 out 104 participants,1.3%). In our cohort (Fig. 2A), the allele was significantly more common (6 out of the 38 participants) compared with both the general gnomAD population (Fisher’s Exact test, P = 1.4*10e−7, OR = 28.17, 95% CI [9.98, 65.45]), the gnomAD African population (Fisher’s Exact test, P = 0.001, OR = 5.66, 95% CI [2, 13.03]) and NIHR Bioresource African controls (Fisher’s Exact test, P = 7.6*10e−4, OR = 13.68, 95% CI [2.41, 143.26]).

The variant, p.Arg185His, results in an amino acid substitution in the linker between D1/S2 and D1/S3, which lies within a voltage sensing domain of Na_v1.7 (Fig. 2B). This residue is highly conserved across all voltage-gated sodium channels. Functional studies have shown that p.Arg185His variant does not alter channel gating properties at room temperature; however, it was associated with enhanced resurgent currents and increases excitability when expressed in dorsal root ganglion neurons.³⁶

Intronic variants and the detection of deletions were included in the gene-level analysis, and none were deemed medically actionable or of interest.

The impact of SCN9A p.Arg185His variant on Na_V1.7 channel function is temperature dependent

Non-freezing cold injury is caused by cold environmental exposure and neuropathic pain is worsened by further cold exposure. We hypothesized that there may be an interaction between Na_V1.7 and an environmental trigger, such that cooling magnifies the effect of p.Arg185Hisvariant on channel function. WT and p.Arg185His Na_V1.7 channels were expressed in combination with β1 and β2 subunits in HEK293T cells and Na_V1.7 mediated currents recorded by whole-cell patch clamp (Fig. 2C–H). Changing the temperature from 20°C to 10°C did not significantly affect the half-activation potential of p.Arg185His (Fig. 2C), but significantly shifted the half-inactivation potential for steady-state fast inactivation of p.Arg185His mutant compared to WT (Fig. 2D). Channel kinetics were slower at 10°C compared with 20°C with similar effects for WT and p.Arg185His (Fig. 2E). We also wanted to assess channel behaviour at higher temperatures (30°C). Increasing the temperature from 20°C to 30°C did not affect I–V curve or the steady-state inactivation of the WT or the p.Arg185His (Fig. 2F and G). Faster inactivation kinetics were observed at 30°C and the change was similar for WT and p.Arg185His (Fig. 2H).

In conclusion, the p.Arg185His mutant exhibited a depolarizing shift of half-inactivation potentials at 10°C but not at 20°C or 30°C, which is a change in channel gating consistent with a gain of function only at lower temperatures.

Gene-wise rare variant association testing results

When comparing European participants against controls, gene-wise, rare variant association testing for Tiers 1–3 genes identified six genes and regulatory regions of four genes with a significant difference in rare variant distribution (Table 7). In total, 177 018 genomic loci were called in Tier 1–3 genes, 132 605 SNPs (93 295 were singletons), 4059 insertions and 8498 deletions (8859 out of the 12 557 were singletons). Variants per loci rate was 0.000046 in priority genes regions. Heterozygosity to Homozygosity ratio was 0.96, 1.35 for SNPs and 0.29 for INDELS. For European participants, gene-wise rare variant association varied according to clinical phenotype (Supplementary Table 2), and 20 genes and 24 regulatory regions reached Bonferroni-adjusted significance for at least 1 phenotype. For African participants, several genes showed similar associations with post-traumatic neuropathy, small fibre neuropathy and non-freezing cold injury (Supplementary Table 2).

For the groups of functionally validated voltage-gated sodium channel variants, a significant association of SCN10A with small fibre neuropathy in Europeans was found. This was driven by two gain-of-function variants, SCN10A(ENST00000449082.2):c.4984G>A, p.Gly1662S ((Allele frequency 0.0039 in Neuropathic Pain and 0.0005 in controls) and ENST00000449082.2:c.1661T>C, p.Leu554Pro (Allele frequency 0 in Neuropathic Pain and 0.00008 in controls) (Table 7). These variants are reported in patients with small fibre neuropathy, cause a gain of function of the Na_V1.8 channel and enhance dorsal root ganglion neuron excitability.^38,39 However, both variants were characterized as likely benign in ClinVar and their allele frequencies in gnomAD suggest that they cannot be the only cause of neuropathic pain. No significant associations were found for SCN9A, nor SCN11A.

For the whole gene set analysis, 146 genes reached Bonferroni corrected significance in Europeans with neuropathic pain when compared to controls (Supplementary Table 3). Four of the genes, KIF1A, KCNQ5, KCNK4 and NOS2 are linked to human neuropathic pain.

The TRPA1 p.Ala172Val variant is associated with episodic pain and gain of channel function

Variants in Tier 2 and 3 genes were filtered (Tables 3–5), using the same approach as for Tier 1 genes to identify those that may be pathogenic. A rare TRPA1 variant, p.Ala172Val (gnomAD allele frequency 0.0003), was identified in a participant who suffers from episodic widespread chronic pain with neuropathic characteristics particularly affecting the trunk. Clinical phenotype was similar to a neuropathic pain syndrome associated with a TRPA1 channelopathy (p.Asn855Ser variant)³²; however, the participant’s pain was not precipitated by physiological stressors. MRI of the brain and spine, nerve conduction studies and skin biopsy of the lower leg were within age- and gender-appropriate reference ranges. The participant’s child, carrying the same variant, was similarly affected by chronic abdominal and pelvic pain with normal investigations. In silico tools, Polyphen and SIFT scores, predicted p.Ala172Val to be deleterious. The variant is moderately conserved across species (Fig. 3A), is situated within the ankyrin repeat domain of the protein (Fig. 3B) and is a non-polar to non-polar amino acid exchange (Grantham score 64). Based on the rarity of the variant, clinical phenotype consistent with the TRPA1 Familial Episodic Pain Syndrome, positive family history and in silico analysis showing the variant to be in the poorly understood N-terminal ankyrin repeats, the TRPA1 p.Ala172Val variant was prioritized for functional studies. These were compared to the TRPA1 variant p.Asn855Ser which is the only variant previously linked to this disorder.

The biophysical properties of WT, p.Ala172Val and p.Asn855Ser hTRPA1 were compared (Supplementary Fig. 3). In HEK293T transfected cells, TRPA1 current–voltage relationship and half-maximal activation for WT, p.Ala172Val and p.Asn855Ser were not statistically different. Under control conditions, WT, p.Ala172Val and p.Asn855Ser current traces showed sustained outward rectification at positive potentials, consistent with an underlying voltage dependence of channel gating. As current density did not differ between WT and the variants, under control conditions, it is unlikely that channel trafficking is affected in p.Ala172Val and p.Asn855Ser channels.

Voltage-dependence activation of hTRPA1 p.Ala172Val and p.Asn855Ser was measured in response to mustard oil (AITC, a TRPA1 electrophilic agonist) and Menthol (non-electrophilic agonist of TRPA1). In response to 25 µM AITC, in the presence of extracellular calcium, p.Ala172Val showed a pronounced linearization of the current–voltage relationship compared to WT in response to a two-voltage-step protocol. A steeper activation curve was observed at positive potentials for both p.Ala172Val and p.Asn855Ser compared with WT (Fig. 3C). Currents were significantly increased at positive potentials for p.Ala172Val (Fig. 3D) and at negative potentials for p.Asn855Ser (Fig. 3E). Analysis of tail currents demonstrated a significant leftward shift of voltage dependence of channel activation in the presence of AITC for both variants (Fig. 3F); however, the slope of the voltage-activation curve did not significantly change. Current–voltage relationships, tested with a voltage-ramp protocol, showed an increase in current densities at positive potentials after the application of 25 µM AITC (Fig. 3G). In summary, these findings show an enhanced response of p.Ala172Val to 25 µM AITC, suggesting a gain-of-function behaviour of this variant, under agonist stimulation.

Application of 100 µM Menthol did not significantly change current density nor voltage sensitivity when applied to p.Ala172Val or p.Asn855Ser in the presence of extracellular calcium (Supplementary Fig. 4). However, in the presence of both extracellular and intracellular calcium, significant changes were observed (Fig. 4A–D). Both p.Ala172Val and p.Asn855Ser channels showed an increase in outward currents and voltage sensitivity in response to 100 µM Menthol.

In summary, the p.Ala172Val variant confers gain-of-function properties on TRPA1 channel in response to the agonists AITC and Menthol. For the latter, the effect was dependent on intracellular calcium.

Discussion

In this study, we applied WGS to a cohort of patients with neuropathic pain disorders. We carried out two analyses. First, we identified medically actionable variants and second, we carried out group-wise association tests at the gene level for rare variants. The diverse phenotypes ranged from congenital insensitivity to pain to painful neuropathy. Clinically relevant findings in genes associated with pain were reported in 12% of participants. The majority of clinically relevant variants were in voltage-gated sodium channels. We made new genotype–phenotype associations, such as the Na_V1.7 p.Arg185His variant which was more frequent in Africans with non-freezing cold injury. We provide new mechanistic insights showing that p.Arg185His interacts with cold, causing a gain of function in the Na_V1.7 gating properties. The gain of function of p.Ala172Val in TRPA1 in response to agonists extends our knowledge of painful TRPA1 channelopathies.

The study of genetic neuropathic pain disorders poses challenges in describing clinical phenotype and assigning pathogenicity to associated variants. We used the gold standard grading system, of the Neuropathic Pain Special Interest Group of IASP.²⁷ Such an approach works well in disorders where there is structural injury to sensory neurons. For example, in painful distal symmetrical polyneuropathy, pain and sensory signs are found in a neuroanatomically plausible distribution (meeting probable criteria), and specialized investigations confirm a lesion of the somatosensory nervous system (meeting definite criteria). However, the majority of Mendelian pain disorders are sensory neuron ion channelopathies, in which pain is episodic, with no structural injury to sensory neurons. For example, in Inherited Erythromelalgia (Na_V1.7) or Familial Episodic Pain Syndrome (TRPA1, Na_V1.9), sensory examination (between pain episodes), neurophysiology and cutaneous innervation are normal. Careful attention to clinical history is, therefore, essential.

We identified medically actionable variants in 12% of the participants. In our cohort, five participants carry pathogenic or likely pathogenic variants. Four pathogenic variants were in voltage-gated sodium channels, which can impact patient care. For example, a pair of sisters with Inherited Erythromelalgia carry the pathogenic SCN9A p.Ile848Thrwith autosomal dominant pattern of inheritance. Only a minority of patients with Inherited Erythromelalgia possess SCN9A mutations.⁴⁰ Identification of a genetic cause for erythromelalgia means family genetic counselling and preferential treatment with non-selective sodium channel blockers,⁴¹ which is not the standard treatment for other causes of neuropathic pain.⁴² The family may also access future treatments such as selective sodium channel blockers.⁴³

A further 20 participants carry VUS in voltage-gated sodium channels that were deemed clinically relevant. Ascribing pathogenicity to ion channel variants is difficult. The majority are relatively common, exhibit subtle channel gain-of-function effects and likely interact with environmental factors. For example, in our cohort the SCN9A(ENST00000409672.1): c.2215A>G, p.Ile739Valvariant was found in four participants with painful neuropathy and described previously as pathogenic.³⁶ However, p.Ile739Valis common in the general population (0.2% allele frequency), while the prevalence of small fibre neuropathy is rare (∼50 per 100 000⁴⁴). It is unlikely that this variant is fully penetrant (thus pathogenic); but may act as a risk factor. Nevertheless, it is important to identify such variants for treatment selection, because Lacosamide (a sodium channel blocker) has shown efficacy in a clinical trial of patients with small fibre neuropathy and rare SCN9A⁴⁵ variants. Voltage clamp analyses can add valuable experimental evidence to the pathogenicity assessment that cannot be replaced by in silico prediction tools, as in silico analysis is not comprehensive and variant’s effect prediction software are likely to underestimate the impact of the gain-of-function effects. However, recent advances in ion channel modelling have improved outcomes.⁴⁶

WGS is now integrated into the UK national health system (see https://www.genomicsengland.co.uk/) and is now considered the standard of care across the UK.⁴⁷ As of 2022 all testing at NHS genomic centres for painful neuropathies and channelopathies is via WGS. The panel application ‘Inherited neuropathies or pain disorder v1.36’ was informed by the NIHR Bioresource (100 000 genomes), our functional studies and includes all the tier 1 pain genes in this paper and TRPA1. Such an approach will identify many more new variants and the findings of our study will be of direct relevance to practitioners assessing patients with neuropathic pain. Understanding how these novel variants relate to clinical phenotype presents significant challenges for clinicians. The combined expertise of our multidisciplinary team meetings was vital in the interpretation of variant pathogenicity. We note that the proportion of cases in which we have found a pathogenic variant that explains the clinical pain phenotype was low (2.4%), although in a higher proportion (a further 10%) the finding was clinically actionable. In many cases, this was VUS in which the relationship to the clinical phenotype could not be causally established, but could still have implications for treatment and so was reported. The low yield in solving cases at a diagnostic level partly represents the fact that neuropathic pain is monogenic in only a minority of cases. In many cases it is likely to arise from gene–environment interaction and/or multiple genes. To try and improve diagnostic yield in the future we could enrich analysis pipelines (for instance to include analysis of repeat expansions in the WGS data), expand co-segregation analysis and also functional analysis. It is not possible to test all novel variants with current patch clamp technology. Automated prediction of pathogenicity and higher throughput functional assays should be a priority, and important next steps in the integration of WGS into healthcare.

We were able to extend genotype–phenotype associations of SCN9A. The SCN9A p.Arg185His variant was more common in study participants with non-freezing cold injury. p.Arg185His is associated with neuropathic pain disorders, such as small fibre neuropathy¹² and painful diabetic neuropathy.¹³ Non-freezing cold injury is a chronic neuropathic pain disorder caused by an acquired sensory neuropathy³⁷ after prolonged exposure to cold and wet environments. Further cold exposure intensifies the neuropathic pain in a similar manner to cold allodynia after platinum-based chemotherapy.⁴⁸ The mechanism for this cold allodynia is unknown. We show that p.Arg185His displays gain-of-function characteristics at lower temperatures. This occurs through a shift of fast inactivation at 10°C, but not at 30°C or room temperature (20°C). Pathogenic SCN9A variants that cause inherited erythromelalgia demonstrate temperature sensitivity at biophysical and clinical levels such that warming intensifies pain and cooling is analgesic.⁴⁹ A link, therefore, exists between the biophysical properties of Na_V1.7 channel, clinical phenotype and temperature sensitivity. An increase of p.Arg185His variant excitability on cooling may contribute to cold allodynia in a sub-group of non-freezing cold injury patients. Another consideration is that increased activity of p.Arg185His at cold temperatures could injure nerves through calcium dependant excitotoxicity.^50,51 In summary, p.Arg185His may contribute to clinical phenotype or increase the risk for cold-induced neuropathic pain, but we cannot conclude that is causative in isolation.

Bevimed is an inference procedure for identifying loci associated with rare hereditary disorders using Bayesian model comparison and was used in an analysis of all cohorts of NIHR BioResource Rare Diseases.¹⁹ Under different modes of inheritance there was strong evidence for 95 genes and 29 binary disease tags (cases versus controls). No diseases from our cohort were among the associations with strong evidence. We used SKAT-O in our cohort. Group-wise association tests can increase the power to detect associations between rare alleles and phenotypes by aggregating variant counts on the gene level that can have an impact on gene function. Neuropathic pain can arise as a complex effect of rare variants with causal effects in the same direction, which the burden test is best powered to detect, or effects with different levels of causality and direction, which the sequence kernel association test is best powered to detect. We used an approach that maximizes test power by considering the combination of both the burden and sequence kernel association tests. In testing for associations for the whole gene set in Europeans with neuropathic pain, 146 genes reached statistical significance. These are gene-wise results showing a significant association of a configuration of rare alleles grouped at the gene level with certain neuropathic pain phenotypes. In most of the associations, the variant-component test was better powered. This is indicated by the low Rho statistics empirically calculated from the data and is suggestive of the presence of variants with different values of causality and opposing direction of effects towards the phenotype, i.e. both protective and deleterious variants present and driving the association. We focused on bi-allelic variants and indels of high quality, with moderate to high impact on protein-coding genes, as we had a moderate sample size and did not want to diminish statistical power after Bonferroni correction. The genes KCNK4, KCNQ5 and NOS2 were significantly associated with neuropathic pain after gene-wise rare variant association testing. These have not been previously linked to human pain at a genetic level, but preclinical models have implicated these genes in pain pathogenesis.^52-54KIF1A, a Tier 2-3 priority gene also reached genome-wide significance. The association with neuropathic pain would require replication in independent cohorts and functional studies.

We undertook an analysis focusing on genes known or likely to be implicated in neuropathic pain, and we, therefore, used more relaxed criteria in filtering rare variants, including modifier, moderate and high-impact variants. We also note that our in silico variant effect predictions are limited by their tendency to underestimate gain-of-function effects. Our analysis showed an association of five Tier 2 genes, including TRPA1, with neuropathic pain. TRPA1 is a calcium-permeable non-selective cation channel that acts as a sensor of noxious external stimuli, such as mustard oil (AITC) and Menthol. The p.Asn855Ser TRPA1 variant is associated with familial episodic pain syndrome characterized by truncal pain triggered by physiological stress or exercise.³² We identified a study participant with a similar clinical phenotype of truncal pain who carried the p.Ala172Val TRPA1 variant. p.Ala172Val is a missense mutation in the fourth domain of the Ankyrin region. Our in vitro electrophysiological studies show that TRPA1 p.Ala172Val variant causes a gain of function in response to agonist stimulation and provides evidence for Ca²⁺-mediated channel activation through the Ankyrin repeat domains of TRPA1 channel.

TRPA1 consists of a large intracellular NH2 and COOH termini, with the NH2 terminus containing an elongated ankyrin repeat domain which is highly conserved. Human TRPA1 consists of 16 ankyrin repeat domains that connect to the transmembrane domains via a linker region. The p.Ala172Val variant did not change the biophysical properties of TRPA1 channel in the naïve state, but did enhance responses to AITC and Menthol. Current density was larger for AITC when compared with Menthol. Due to differences in their electrophilic nature, agonists vary in their ability to covalently (AITC) or noncovalently (Menthol) modify the channel upon binding. This difference in channel binding might underlie different agonist responses. The p.Asn855Ser and p.Ala172Val variants increase activation in response to agonists, although with distinct impacts on the biophysical properties of the channel. Intrinsic differences due to the positions of the variants within the channel (p.Asn855Ser is in the transmembrane domain S4), might underlie the differences observed. Our data show that the ankyrin repeat domain is important in TRPA1 channel activation and gating functions. The p.Ala172Val variant is reported at a heterozygous frequency of 0.0003 in gnomAD and individuals over 70 years of age are reportedly healthy so it is unlikely to be fully penetrant. It is more likely to act as a risk factor and contribute to the expression of neuropathic pain by altering the functional properties of nociceptive afferents; however, the drivers for channel activation (environmental versus endogenous ligands) are not known.

In our study, we identified clinically relevant variants in 12% of the participants, with an impact on clinical care and treatment. The majority of these variants are located in ion channels which are enriched in nociceptors and environmental triggers (such as cold in the case of non-freezing cold injury) may interact and enhance the gain-of-function impact of such variants.

Supplementary material

Supplementary material is available at Brain Communications online.

Acknowledgements

We thank NIHR BioResource volunteers for their participation, and gratefully acknowledge NIHR BioResource centres, NHS Trusts and staff for their contribution. We thank the National Institute for Health and Care Research, NHS Blood and Transplant, and Health Data Research UK as part of the Digital Innovation Hub Programme. The views expressed are those of the authors and not necessarily those of the NHS, the NIHR or the Department of Health and Social Care.

Funding

This research was funded in whole, or in part, by the Wellcome Trust (Grant numbers: 109915/Z/15/Z, 083259, 202747/Z/16/Z). For the purpose of open access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.

D.L.H.B., A.C.T. and A.S.C.R. are members of the DOLORisk consortium funded by the European Commission Horizon 2020 (ID633491), which received funding from European Union’s Seventh Framework Program (FP7/2007-2013). D.L.H.B. and A.C.T. are members of the International Diabetic Neuropathy Consortium (IDNC) research programme, which is supported by a Novo Nordisk Foundation Challenge Programme grant (Grant number NNF14OC0011633). D.L.H.B. is a Wellcome Investigator (202747/Z/16/Z and 223149/Z/21/Z). This work was supported by the Wellcome Trust through a Strategic Award to the London Pain Consortium (ref. no. 083259). D.L.H.B. and A.C.T. are members of the PAINSTORM consortium funded by UK Research and Innovation (UKRI) and Versus Arthritis. A.C.T. is supported by Academy of Medical Sciences Starter Grant SGL022\1086 and is an Honorary Senior Research Fellow and Carnegie-Wits Diaspora Fellow at the Brain Function Research Group, School of Physiology, Faculty of Health Sciences, University of the Witwatersrand, Johannesburg, South Africa. R.H. is a Wellcome Trust Investigator (109915/Z/15/Z), who receives support from the Medical Research Council (UK) (MR/N025431/1 and MR/V009346/1), the European Research Council (309548), the Newton Fund (UK/Turkey, MR/N027302/1), the Addenbrookes Charitable Trust (G100142), the Evelyn Trust, the Stoneygate Trust, the Lily Foundation and Medical Research Council (MRC) strategic award to establish an International Centre for Genomic Medicine in Neuromuscular Diseases (ICGNMD, MR/S005021/1). M.M.R. acknowledges funding from the Medical Research Council (MRC MR/S005021/1), the National Institutes of Neurological Diseases and Stroke and Office of Rare Diseases (U54NS065712 and 1UOINS109403-01 and R21TROO3034), Muscular Dystrophy Association (MDA510281) and the Charcot Marie Tooth Association (CMTA). This research was supported by the National Institute for Health Research (NIHR) Cambridge Biomedical Research Centre (BRC-1215-20014). D.L.B and A.C.T were funded by the National Institute for Health Research (NIHR) Oxford Biomedical Research Centre (BRC).

Competing interests

D.L.H.B. has acted as a consultant on behalf of Oxford Innovation for Amgen, Bristows, LatigoBio, GSK, Ionis, Lilly, Olipass, Orion, Regeneron and Theranexus on behalf of Oxford University Innovation over the last 2 years. He has received research funding from Lilly. He has received an industrial partnership grant from the BBSRC and AstraZeneca.

ASCR declares the following interests occurring in the last 24 months: undertakes consultancy and advisory board work for Imperial College Consultants that included remunerated work for Abide, Confo, Vertex, Pharmanovo, Lateral, Novartis, Mundipharma, Orion, Shanghai SIMR Biotech Asahi Kasei, Toray & Theranexis; owner of share options in Spinifex Pharmaceuticals from which personal benefit accrued upon the acquisition of Spinifex by Novartis in July 2015 (final payment was made in 2019); named as an inventor on patents: Rice A.S.C., Vandevoorde S. and Lambert D.M Methods using N-(2-propenyl)hexadecanamide and related amides to relieve pain. WO 2005/079771 and Okuse K. et al. Methods of treating pain by inhibition of vgf activity EP13702262.0/WO2013 110945; Member Joint Committee on Vaccine and Immunisation—varicella sub-committee; Analgesic Clinical Trial Translation: Innovations, Opportunities, and Networks (ACTTION) steering committee member; Medicines and Healthcare Products Regulatory Agency (MHRA), Commission on Human Medicines—Neurology, Pain & Psychiatry Expert Advisory Group.

Data availability

The genotype and phenotype data can be accessed by application to the NIHR BioResource Data Access Committee at [email protected] or by application to Genomics England Limited following the procedure outlined at https://www.genomicsengland.co.uk/research.

References

Abbreviations

 
• AITC =

  Allyl isothiocyanate

 
• HGNC =

  HUGO Gene Nomenclature Committee

 
• HSAN =

  Hereditary and Autonomic Sensory Neuropathy

 
• MAF =

  Minor allele frequency

 
• NFCI =

  Non-freezing cold injury

 
• NIHR =

  National Institute for Health and Care Research

 
• NOS =

  Not otherwise specified

 
• NPD =

  neuropathic pain disorders

 
• SNP =

  single nucleotide polymorphism

 
• VUS =

  variant of uncertain significance

 
• WT =

  Wild type

Author notes