Introduction
Mevalonate kinase deficiency (MKD; 610377) is an autosomal monogenic recessively-inherited disease, caused by mutations in the MVK gene (12q24.11) coding for mevalonate kinase (MK). MK is a key enzyme of the mevalonate pathway, which is essential for the biosynthesis of isoprenoids and the decrease of which is considered to lead to an overproduction of the specific marker of disease, interleukin (IL)-1β (1–3). These mutations lead to a shortage of intermediate compounds and final products of the metabolic route of cholesterol (4–6).
The residual activity of MK defines different degrees of MKD severity, ranging between an auto-inflammatory phenotype (hyper IgD syndrome/HIDS; OMIM #260920) and severe clinical presentation (mevalonic aciduria/MA; OMIM #610377) (6). HIDS is characterized by fever, which re-occurs every 2–8 weeks, variably accompanied by malaise, headache, diarrhea, abdominal pain, vomiting, skin rash, arthralgia, arthritis, tender lymphadenopathy, anemia normocytic, hepatosplenomegaly, and as oral and genital ulcers. Urinary mevalonic acid (UMA) is increased predominantly during inflammatory episodes, whereas increased levels of IL-1β, IL-1α, IL-6, IL-10, IL-18, tumor necrosis factor (TNF)-α and IFN-γ cytokines are also observed between fever episodes (7,8). Laboratory analyses show increased neutrophil counts and high acute phase reactants, often associated with persistently high serum IgD and IgA (9).
Patients with MA exhibit the above-mentioned symptoms associated with dysmorphic features, cataracts, uveitis, neurological impairments, and failure to thrive. UMA excretion is also high between attacks (10). There have been >70 different disease-causing mutations identified in the MVK gene. The majority of these exhibit poor genotype-phenotype correlation with respect to clinical presentation and biochemical derangements (11). Notably, it has been found that >95% of patients with HIDS are compound heterozygous for the V377I MVK allele, whereas a second mutant allele, I268T, is specific to patients with MA (12). Other mutations have been described in patients with HIDS and MA without a reliable genotype/phenotype correlation (13), although patients carrying the same mutation often exhibit substantial variability in symptoms and respond differently to treatment (14).
Of note, previous studies have described a correlation and/or an overlap between the MKD phenotype and/or the MVK genotype with other diseases, including inflammatory bowel diseases (15,16), possibly due to shared genetic background (17), and retinitis pigmentosa, adding novel features to the wide range of the already known MKD phenotypes (18,19).
Whilst several reports suggest that the shortage of isoprenoid intermediates, including geranylgeraniol, is involved in the autoinflammatory phenotype of MKD, further data has led to a focus of attention on the end products of this pathway, including 25-deoxycholesterol (CH25) (20–22). Of note, 25-hydroxylcholesterol can act by modulating the mevalonate pathway itself by suppressing sterol regulatory element-binding proteins (SREBPs), thus highlighting a complex feedback network, which may be relevant in the modulation of MKD phenotypes (23,24). The present study hypothesized that MKD may require consideration a polygenic, rather than a monogenic, disease. The overlap and the wide range of clinical signs shared by MKD and other auto-inflammatory diseases can cause delayed or incorrect diagnosis. In order to verify this hypothesis, the presents study performed whole exome sequencing (WES) analysis on DNA samples obtained from patients with mild or severe MKD phenotypes to determine any more or less pathogenetic variants, which may be involved in determining the phenotype of the patient.
Materials and methods
Study design
Initially, the present identified the small nucleotide variants (SNVs)/small insertion/deletions (INDELs), which had a minor allele frequency (MAF) <0.05, as reported in the NHLBI Exome Sequencing Project (ESP) Exome Variant Server (http://evs.gs.washington.edu/EVS/) database and referred to the general population. In order to be considered, the variants also had to be harbored by genes belonging to the cholesterol biosynthetic pathway (Table I). The hypothesis underlying this approach was that the variants harbored by these genes modulate overall cholesterol biosynthetic activity, possibly affecting the overall availability of cholesterol and/or its biosynthetic intermediates, thus contributing to the modification of MKD clinical phenotypes.
The present study also investigated an alternative analytical strategy, selecting only the variants that are poorly represented in the general population and, at the same time, predicted as potentially pathogenic. For this purpose, the SNVs/INDELs were selected, according the following inclusion criteria: (a) variants with an MAF <0.03 in the general population, as reported in the ESP database, (b) SNVs leading to a non synonymous amino acid substitution, (c) SNVs/INDELs carried in homozygous state, (d) SNVs/INDELs predicted as pathogenic by the following in silico algorithms: Polyphen-2 (25), Mutation Taster (26) and likelihood relation test (27), according to scores recorded in the dbNSFP v2.0 database (http://sites.google.com/site/jpopgen/dbNSFP); (e) SNVs/INDELs considered phylogenetically conserved, based on genomic evolutionary rate profiling (GERP)++ scores reported in the dbNSFP v2.0 database. The first analytical approach examined the hypothesis that a number non-rare variants are involved in modifying the MKD clinical phenotype, acting along the same biochemical pathway, whereas the second was intended to identify rare and pathogenic variants, carried in a homozygous state and possibly associated to certain, more severe, MKD clinical phenotypes.
The two strategies were performed taking into consideration the phylogenetic nucleotide evolutionary conservation, based on PhyloP (28) and GERP++ scores, as reported in dbNSFP v2.0; and only substitutions of conserved nucleotides were considered for their potential pathogenic role. This type of alteration is a predictor of deleteriousness, being a variation that reduces organism fitness, which is a property closely associated with molecular pathogenicity (29).
DNA extraction
Genomic DNA (gDNA) was extracted from 1–2 ml EDTA-anticoagulated blood from the proband and their parents using an EZ1 DNA Blood kit (Qiagen, Milan, Italy), according to the manufacturer's protocol.
Whole exome sequencing and bioinformatics analysis
The technical and scientific review board of the Institute for Maternal and Child Health-IRCCS 'Burlo Garofolo' (Trieste, Italy; no.185/08; 19/08/2008) approved the present study. For a child to be eligible, informed consent had to be obtained from parents or caregivers. Furthermore, patients with MKD of any age were excluded from the study if they had an acute or chronic infectious disease, any clinically significant disorder, or if they were currently on any medication with known effects on immunological factors, including corticosteroids. Blood was collected by venepuncture from five patients with mild or severe MKD phenotypes, determined by a physician using a visual analogue scale score (Table II). For each patient, a specific control of the identical sex and the age was used. Each control sample was collected, preserved and analyzed in an identical method of the patient to precisely identify phenotypic differences.
The samples obtained from the patients diagnosed with MKD were analyZed using whole exome sequencing (WES). Starting with 3 µg gDNA, collected from each subject, a TruSeq™ Exome Enrichment 62 Mb kit (Illumina, Inc., San Diego, CA, USA) was used to capture the whole exome; of which the target size was 62 Mb, comprising the overall genes coding sequence, the 5′-untranslated regions (UTRs) and 3′-UTRs. A fragment exome library was constructed, according to the manufacturer's protocol, and a 100 bp paired-end sequence was analyzed using the Illumina HiSeq 1000 platform (Illumina, Inc.). Raw sequencing data were collected as unmapped reads in fast Q format. CLC Genomics Workbench ver. 6.5 software was used to assess the quality of reads, to map reads back to the human reference genome, hg19, to calculate the overall coverage, to perform local realignment and base quality recalibration, and to identify SNVs and INDELs, all of which were collected into standardized Variant Call Format version 4.1 (30). The SNVs/INDELs were annotated using ANNOVAR software (31) referring to the following public databases: refGene, NCBI dbSNP build137 (http://www.ncbi.nlm.nih.gov/SNP/), 1000 Genomes Project (http://www.1000genomes.org/), NHLBI ESP Exome Variant Server (http://evs.gs.washington.edu/EVS/), dbNSFP v2.0 (27) and NCBI ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/).
Polymerase chain reaction (PCR) and sanger analysis
PCR amplification was performed for the specific gene coding sequence of the MVK and RAB40AL genes in 15 µl total volume. Each reaction contained 50 ng/µl genomic DNA and KAPA 2 G Fast Hot Start Readymix (RESNOVA, Rome, Italy). PCR amplification was performed using a common annealing temperature in a touchdown thermocycler with a two-step cycle: Initial denaturation at 96°C for 3 min, followed by a gradiant (0.5°C reduction per cycle) of 10 cycles at 96°C for 15 sec, 63°C for 10 sec and 72°C for 1 second. The second step was 28 cycles of 96°C for 15 sec, 53°C for 15 sec and 72°C for 1 sec. The PCR products were purified with 2 µl ExoSAP (USB Corporation, Cleveland, OH, USA) by incubation at 37°C for 20 min and 85°C for 10 min. All thermal cycling, PCR amplifications, ExoSAP purifications and sequencing reactions were performed in a 2720 Thermal Cycler (Applied Biosystems; Thermo Fisher Scientific, Inc., Waltham, MA, USA). The amplification products were directly sequenced by the Sanger method using an ABI PRISM 3130XL automated DNA sequencer (Applied Biosystems; Thermo Fisher Scientific, Inc.). The sequences were analyzed using Seqman II software (DNASar I Lasergene, 7.0; DNAStar, Inc., Madison, WI, USA). The PCR amplification and Sanger sequencing primers used to confirm the variations identified in the PEX11 G gene were as follows: PEX11GcDNA, forward 5′-TGA AAC TGA GAC AGA GGCTG-3′ and reverse 5′-AGT GTC AGG GGG TAG TGG-3; and Pex11gcDNA, forward CCT GTG GAC AAT GCT GAAG and Pex11gcDNA, reverse TCA TCA AGG GCT GTC TGC.
The PCR amplification and Sanger sequencing primers used to confirm the variations identified in the uracil-DNA glycosylase (UNG) gene were as follows: UNG exon2, forward 5′-CTG TCC GCT TTT GCT GGG-3′ and reverse 5′-CCG GCT ACA CTA ACA AGAC-3′.
Urinary mevalonic acid measurement
Mevalonic acid levels were determined in the patients with MKD using standard procedures, as described by Shoemaker et al (32). The urine samples were collected over 24 h and stored at 4°C until analyzed.
Urinary mevalonic acid concentrations were determined using gas chromatography/mass spectrometry. Spectra were obtained using a Hewlett Packard gas chromatograph 6890 system (Hewlett Packard, Palo Alto, CA, USA), equipped with a Hewlett Packard 5973 quadrupole, operating in electron-impact mode at 70 eV, as described previously (32). UMA is considered to be within a normal range between 34 and 323 µg/ml (internal laboratory reference).
B lymphocyte phenotype
The analysis of B lymphocyte phenotype was performed using the whole peripheral blood samples from the patients and controls. Briefly, 100 µl of blood was washed twice in phosphate-buffered saline (PBS), and then stained with anti-IgD/IgM fluorescein isothio-cyanate (2.5 µl/test; cat. no. 348206/314506; BioLegend, Inc., San Diego, CA, USA), anti-CD27 phycoerythrin (PE; 5 µl/test; cat. no. 130093185), anti-CD19 PE-Vio770 (5 µl/test; cat. no. 130096641); and anti-CD45 (5 µl/test; cat. no. 130092880), all from Miltenyi Biotec GmbH (Bergisch Gladbach, Germany). Data were acquired using a Cyan ADP cytometer (Beckman Coulter, Brea, CA, USA), and analysis was performed using FlowJo software v7.6 (Tree Star, Inc., Ashland, OR, USA). B-switched memory cells were identified as CD45/CD19/CD27 positive and IgM/IgM negative.
Results
A total of ~14 Gb of sequence data per sample were produced using WES analysis, corresponding to 53X overall median coverage, among all the samples. The overall target percentage covered at least 20X, and was 80.7, 71, 73.6, 79, 75.3% of the whole exome, in patients 1, 2, 3, 4 and 5, respectively. On average, 59.404 variants per sample were identified.
As shown using Sanger Sequencing, the patients with MKD in the present study carried at least one mutation in the MVK gene. Data are present, together with the patients' mevalonic acid levels, in Table III.
By selecting the SNVs/INDELs with an MAF <0.05, the present study identified 10 exonic SNVs, listed in Table IV, six of which were harbored by MVK, as already known, with the remaining four variants identified in the SQLE, IDI1 and NSDHL genes, even when none of them were in a homozygous or a compound heterozygous state.
Subsequently, the SNVs/INDELs were selected, according to the second filtering criteria, which can define SNPs that are hard to confirm. The results of this filtering identified four variants in patient 2 (P2) and two variants in patient 5 (P5), as shown in Table V.
Variants in P2
The analyses revealed an SNV in the PEX11γ gene (19p13.2), carried in a homozygous state by patient P2. A missense variation, rs11668511 (NM_080662 c.C646T; p.L216F), was found in exon 5 of the PEX11γ gene. Several in silico tools, including Polyphen-2, MutationTaster and LTR, predicted this variant to be pathogenic or disease causing (Table V).
The present study also found the following variants in different genes: rs146931399 (NM_145246.4 c.G241A; p.V81I) on the FRA10AC1 gene, rs145520946 (NM_031302.3 c.A817G; p.M273 V) on the GLT8D2 gene and, finally, rs202215334 (NM_001080489 c.G304A:p.D102 N) on the GLOD5 gene. The functions of the proteins encoded by these genes remain to be elucidated and limited previous data exists. Therefore, it is difficult to correlate these variants with the clinical phenotype of P2.
Variants in P5
The present study detected a variant in the UNG gene (12q23-q24.1), which was carried in a homozygous state in patient P5. In this patient, the missense variation, rs151095402 (NM_080911.2 c.262C>T; p.R88C), in exon 2, was found in homozygosis for the first time, to the best of our knowledge.. Despite being rare (MAF, ~0,0012), this missense mutation had no effect on the IgG, IgA or IgE concentrations. In addition, a normal percentage of B-switched memory lymphocytes was observed (Fig. 1).
Discussion
Several studies have documented the phenotypic heterogeneity of patients carrying mutations in the MVK gene (15–19).
In the patients recruited in the present study, a correlation was observed between MVK genotype and UMA urine levels, and, of note, p.V132I has already been (33) described as being associated with lower UMA levels. However, the present study found a poor correlation between the' clinical phenotypes of the patients and MVK mutations, even taking into account the mutations, which change amino acids located close to each other in the protein sequence, for example, p.V377I and p.G336S, and that were considered to affect the same MK protein domain. Of considerable clinical importance is the p.G336S mutation which, unlike other mutations, determines inhibition of the metabolic pathway, causing high UMA levels, which also occur outside of fever attacks (5).
With the aim of identifying novel potential modifier genes/variations associated with MKD phenotype variability, the present study analyzed the entire exome of the five patients with MKD, hypothesizing that MKD may be a polygenic disease. Considering the heterogeneity of the clinical phenotypes, the present study investigated the possible presence of modifier genes with potential to affect the phenotype of the patients, in any way.
The variants obtained from a first analysis, which filtered variants with an MAF <0.05 did not explain the wide range of MKD phenotypes observed in the patients
When the analysis was performed with the second filtering criteria, two putative genes were identified, which potentially explained the clinical phenotype of two MKD patients (P2 and P5).
Patient P2 was identified with a mutation in the PEX11γ gene, which is a member of the PEX11 family, whose predominant function is associated with the tabulation, enlargement and clustering of peroxisomes, thus being important in human metabolism (34). Peroxisomal gene defects are also known to be associated with severe disorders (35–38).
In mammals, there are three PEX11-associated genes, PEX11-α, PEX11-β and PEX11-γ. Several previous studies analyzing expression have been performed on PEX11-α, and a few studies using animal models have been performed examining the β proteins (35,39). The γ protein has been reported to interact with α and β, and, by being present at the crossroads of PEX11, activates the peroxisome proliferation pathways, possibly being involved in homotypic interactions. In addition, PEX11-γ overexpression has been suggested to have the ability to induce early formation of juxtaposed elongated peroxisomes, suggesting that the γ protein either acts upstream of PEX11-α and PEX11-β, or that is their limiting factor. PEX11-γ appears to recruit the other PEX11 proteins to facilitate aggregation, and is required for the elongation of the peroxisome membrane (40).
Previously, a case of a patient carrying a homozygous nonsense mutation in the PEX11-β gene has been described. Using immunofluorescence microscopical analyses in over-expression experiments at 37°C/40°C, it was suggested that the occurrence of the enlarged catalase-containing peroxisomes is associated with the PEX11-γ gene (38).
The PEX11γ gene and the c.C646T; p.L216F variant have received limited attention in previous studies. The present study hypothesized that this variant, rarely detected in the general population (MAF, ~0.024), affected and modified the phenotype of P2. Patient P2 presented with elevated UMA values due to mutations in the MVK gene. However, in addition to showing all the typical signs of a MKD phenotype, the patient presented with visual blurring, not associated with cataracts (37).
A mutation in the UNG gene encoding an important DNA repair enzyme, was also found in P2.
The deletions in the UNG gene (c.391_393delC, c.426_429delAT c.568_571delTA) and the missense mutations (c.752T>C p.F251S) have been previously described to be associated with hyper IgM immunodeficiency type 4 (40,41). Hyper-IgM syndrome is characterized by normal or increased levels of IgM, and is associated with low or absent serum levels of IgG, IgA and IgE, indicating a defect in the class-switch recombination process (42). In 2003, Imai et al (43) reported on three patients with HIGM5, in which deleterious mutations within the catalytic domain of the UNG protein were identified. Functional investigations of immunoprecipitation and immunoblot assays confirmed the absence of expression in the three patients, supporting the hypothesis of a correlation between the identified mutation and protein instability.
In the present study, the p.R88C missense mutation in the UNG gene identified in P5, was previously described by Torseth et al (41), but only in a heterozygotic state. Despite the fact that this mutation is rare, its presence in homozygosis in the patient in the present study had no effect on the levels of IgG, IgA and IgE (Fig. 1). In addition, the patient was found to have a normal percentage of B-switched memory lymphocytes, indicating normal immunity, with regards to the B-cell compartment. Furthermore, patient four (P4) who was not a carrier of mutations in the UNG gene, had high levels of IgA, as previously reported described in several patients with MKD (44,45) but with normal levels of IgG and IgM. Taken together, these results suggested that investigation of the polymorphism found on the UNG gene requires caution, as it does not correlate directly with any specific disease phenotype.
The primary aim of the present study was to identify novel potential modifier genes for MKD disease. A number of preliminary results have shown that the clinical profile of heterogeneous phenotypes of patients with MKD may be associated with novel genes involved in modulating the MKD clinical picture (15,18).
In the present study, following WES in five patients with MKD, it was found that variations in genes encoding proteins of the cholesterol pathway were not associated with the modulation of MKD clinical phenotypes.
A missense variation, namely the c.C646T; p.L216F (NM_080662) in exon 5 of the PEX11γ gene was observed in homozygosis in P2, possibly correlating with visual blurring. A rare UNG gene variant, namely (c.C262T; p.R88C), was detected in homozygosis in P5, however this did not correlate with a specific clinical phenotype. A number of other variants were found in the five MKD patients analyzed, however, no correlations with phenotype were observed. The lack of a direct correlation between genetic variations and phenotype suggest a possible role of post-transcriptional mechanisms, which may affect protein expression or function.
In the patients examined in the present study, no mutation was identified, but a a synonymous variation on NM_003956:exon1:c.C657T:p.N219 N in the gene coding for the CH25H gene was identified in the P2, P4 and P5 patients.
The patients in the present study were also analyzed using direct sequencing to verify the presence of the intronic polymorphism, identified by Moura et al (46). The intronic NM_001510.3:c.89-32007 A>G polymorphism (rs1450500) of the human glutamate receptor δ-2 (GRID-2) gene was not identified in any of the patient. Therefore, it was concluded that, at least in these patients, and in patients with more severe features, GRID2 does not appear to be significantly involved.
Taking into account the primary limitation of the present study lacking patients with MKD exhibiting different clinical phenotypes, but sharing the same MVK mutations, the use of WES, despite being an attractive approach, only provided evidence of an association between one genetic variant, in the PEX11γ gene, and one clinical characteristic, visual blurring, in one of five patients. Therefore, further analysis required, using NGS approaches on a larger sample size of patients with MKD sharing the same MVK mutations and, ideally, exhibiting extreme clinical phenotypes, in order to identify genes and variants that correlate with clinical features. Considering the fact that MKD is a rare, orphan disease, this approach may be possible providing there are a sufficient number of patients with the above-mentioned characteristics. The identification of modifier genes specific for HIDS and MA may assist in the diagnoses of these two forms of the same disease at an earlier stage.