Biallelic variants in NSUN6 cause an autosomal recessive neurodevelopmental disorder

,


Introduction
Intellectual disability (ID) and developmental delay (DD) are neurodevelopmental disorders characterized by highly heterogeneous phenotypes that affect 1% to 3% of the general population. 1,2The introduction of next-generation sequencing technology revealed the association of many single genetic defects, implicating >1700 genes in ID. [2][3][4] Although autosomal dominant de novo variants are the predominant cause of ID in outbred populations, 5,6 autosomal recessive gene defects are the leading genetic cause of ID in countries with frequent parental consanguinity. 7][10][11][12] Consanguineous families thus provide a unique opportunity to identify novel recessive causative genes. 13NA modifications have emerged as an important posttranscriptional regulation mechanism, modulating gene expression of different biological processes. 14Over the years, >170 RNA modifications have been identified in both coding and noncoding RNAs. 15,160][31][32] Different levels of m 5 C in mRNA have been observed among mammalian tissues. 19,33Most of the m 5 C RNA modifications are driven by RNA methyltransferase of the NOL1/NOP2/SUN domain (NSUN) family, a family counting 7 members in human (NSUN1-7). 34Nsun transcripts are developmentally regulated during mouse embryogenesis and show tissue-specific enrichment. 357][38][39][40][41][42][43][44][45] NSUN5 (aka, WBSCR20A) and its 2 truncated paralogs NSUN5P1 and NSUN5P2 (aka, WBSCR20B and WBSCR20C, respectively) map to the region commonly deleted in Williams-Beuren syndrome or its flanking sequences and may contribute to certain of its phenotypical features. 468][49][50] Here, we report on 3 consanguineous unrelated individuals with DD and ID, carrying biallelic deleterious variants in NSUN6.The homozygous deletion of the Drosophila ortholog of NSUN6 led to locomotion and learning impairment.Overall, our data provide evidence that biallelic pathogenic variants in NSUN6 cause one form of ARID and establish another link between RNA modification and cognition.

Variant identification and protein modeling
The exome of affected individual 1 was captured with biotinylated probes (Twist Human Core Exome EF Multiplex Complete kit), following the manufacturer recommendations and paired-end (2 × 75 bp) reads on an Illumina NextSeq 500 sequencer.Sequences were analyzed in accordance with the SeqOne software procedures based on good practice from the Broad Institute (Internal Pipeline: BWA-MEM, GATK v3.6; SeqOne Pipeline: v 1.2, 2018).Variants were filtered according to their level of coverage (>10-fold), their frequency in Genome Aggregation Database (gnomAD) (<1%), their quality, and the predicted impact on the protein.
The exome from the affected individuals 2 and 3 was captured using the xGen Exome Research Panel v2 (Integrated DNA Technologies) and sequenced using the Illumina HiSeq4000 platform according to the manufacturer's protocols.The overall mean-depth base coverage was 145and 138-fold and 98% and 97% of the targeted region was covered at least 20-fold, for individuals 2 and 3, respectively.Read mapping and variant calling were performed as described. 51,52Variants not passing the quality and not predicted to affect the protein were filtered out.Homozygous and hemizygous variants with a minor allele frequency (MAF) < 1% in the general population (1000 Genomes, Exome Variant Server, and gnomAD) were retained (Supplemental Table 1).
Nonsynonymous variants were modeled using the PDB entry 5WWR. 53Chains A (NSUN6) and C (tRNA) were used to highlight the position of the variants with Swiss-PdbViewer. 54

Methylation assay
Histidine-tagged NSUN6 variants and human tRNA transcripts were prepared and purified as previously described. 53,55ethylation assays and enzymatic kinetic measurement methods have been previously described in detail. 53,55Briefly, a reaction mixture consisting of 200 μM 3 H-S-adenosyl-Lmethionine (SAM), 50 mM Tris-HCl, pH 7.0, 100 mM NaCl, 10 mM MgCl 2 , 100 μg/mL bovine serum albumin, 5 mM DTT, and 5 μM tRNAs (GCA) was initiated by the addition of 500 nM enzyme.At time intervals ranging between 2 and 8 minutes, 5 μL aliquots were removed to glass fiber filter disks and soaked in 5% trichloroacetic acid to precipitate the labeled methylated tRNA.After washing, the amount of radioactive [methyl-3 H] tRNA on each disk was measured in a Beckman Ls6500 scintillation counting apparatus.
RNA mass spectrometry analysis tRNA transcripts (200 ng) reacted with NSUN6-WT or NSUN6-D323N were hydrolyzed with 0.2 μL benzonase, 0.25 μL phosphodiesterase I, and 0.25 μL bacterial alkaline phosphatase in a 20 μL solution including 4 mM NH4OAc at 37 • C overnight.After complete hydrolysis, the products were dissolved in acetonitrile and then applied to ultraperformance liquid chromatography-mass spectrometry/ mass spectrometry (UPLC-MS/MS).The nucleosides were separated on a Hilic column (Atlantis HILIC Silica 3 μM, 2.1 × 150 mM column) and then detected by a triple quadrupole mass spectrometer (AB Sciex QTRAP 6500+) in the positive ion multiple reaction monitoring mode.The nucleosides were quantified using the nucleoside-to-base ion mass transitions of m 5 C (Q1/Q3 = 258.1/126.1)and A (Q1/Q3 = 268.1/136.2).Quantification was performed by comparison with the standard curve obtained from pure nucleoside standards running in the same batch.Mobile phase A contains 50% acetonitrile and 50% deionized water.Mobile phase B contains 95% acetonitrile and 5% deionized water.

Isothermal titration calorimetry assay for SAM binding
Isothermal titration calorimetry (ITC) measurements were performed at 25 • C, using MicroCal PEAQ-ITC (Malvern Panalytical) or isothermal titration microcalorimetry ITC-200 (Malvern Instruments).Experiments were performed by titration of 20 injections of 2 μL 1 mM SAM into the sample cell containing around 50 μM purified NSUN6-D323N mutant protein solution.The SAM and corresponding protein were held in the same buffer (50 mM Tris-HCl (pH 7.0) and 100 mM NaCl).Titration of SAM to the same buffer was used as control.Binding isotherms were fitted by nonlinear regression using MicroCal PEAQ-ITC analysis software or Origin Software version 7.0 (Micro-Cal).The ITC data were fitted to a one-site binding model using the 2 software as described.

Behavior assays
For the negative geotaxis assay, 5-day-old male flies were placed 1 day before the experiment in groups of 10 individuals into a fresh food vial.On the day of the experiment, these flies were transferred without anesthesia into measuring cylinders and allowed to adapt to the new environment for 15 minutes.The flies were gently tapped to the bottom of the cylinder, and the number of flies climbing above a 10-cm mark was measured.The assay was repeated for each group 6 times.Four individual groups were used per genotype.
Four days before the courtship learning experiment, male test subjects were collected directly after eclosure and kept isolated in individual vials until the day of experiment.In addition, female Canton-S (Wild-type) virgins were collected, and a standardized premated female culture was established by combining 20 wild-type virgins with 10 wildtype males.On the day of experiment, a premated female was placed to each isolated male, and the cap of the tube was gently pushed into the tube to increase the contact between the flies (training phase).Sixteen naive male flies per genotype were left isolated.After 1.5 hours, the naive and trained male flies were directly transferred into courtship chambers together with a wild-type female virgin, and the flies were recorded for 10 minutes.Sixteen female-male pairs were recorded in parallel, whereas naive and trained male flies were included in each recording to avoid differences because of the circadian time.The time the male spends courting was quantified by observing distinct courtship behavior as orienting toward and following the female, tapping, extending and vibrating the wing, licking, and (attempting) copulation.The videos were randomized to allow a blinded scoring of the male courtship behavior.The courtship index was calculated as Ci = (time spent courting)/(total amount of time) and the learning index as learning index = (Ci(naive)-Ci(trained))/(Ci(naive)).

Patients with biallelic variant in NSUN6
We identified 3 unrelated individuals with biallelic deleterious variants in NSUN6 and presenting with DD/ID through exome sequencing.
Individual 1 is a girl born from first cousin parents from Yemen.She is presenting with severe growth delay, ID, global DD, motor delay, microcephaly, and seizures (Table 1).Specifically, her pregnancy was complicated with oligohydramnios, and she was delivered prematurely at 34 weeks of gestation with a birth weight of 1.4 kg and a head circumference of 30 cm.She was described as dysmorphic and with mild low-set ears and moderate jaundice at birth.A computed tomography scan at 6 days of age reported mild cerebellar atrophy and lissencephaly.A follow-up magnetic resonance imaging at 2 months of age confirmed hypoplastic cerebellar vermis.Sagittal T1-weighted midline image showed prominent cerebrospinal fluid space at posterior fossa.There was incomplete lobulation of the vermis, with absence of the prepyramidal fissure.The vermis was rotated superiorly with widening of the tegmento-vermian angle.Ventricular system and basal cisterns were normal.The signal intensity of the parenchyma was normal, and there was no cortical dysplasia or evidence of subependymal heterotopia.Corpus callosum was normal.No lissencephaly findings were reported.At 14 months of age, her weight, height, and head circumference were below the mean for her age (3 kg, 56 cm, and 34 cm, respectively).She was not walking, standing, seating, or rolling over in bed.She was unable to pronounce any words or syllables.Her face did not look particularly dysmorphic except for her notorious microcephaly.She has 4 unaffected siblings and 1 affected brother who died at 4 years of age with the same symptoms and blindness (Figure 1A).We identified a homozygous frameshift variant by exome sequencing: NM_182543.2:c.25_26del;NP_872349.1:p.(Leu9Glufs*3).Segregation analysis could not be tested.This variant is observed in 2 heterozygous individuals in gnomAD (MAF = 0.000009188) and in the REGENERON Genetics Center database (MAF = 0.000004).It is not present in TopMed nor in the Iranome database.This homozygous frameshift variant maps to the first exon of NSUN6, suggesting that the mutated transcript is likely eliminated by nonsense-mediated decay, leading to the absence of the enzyme.
Affected individual 2 is an 11-year-old boy presenting with moderate ID, speech impairment, global DD, motor delay, and behavioral anomalies (aggressiveness and sleeping difficulties).Height, weight, and occipital frontal circumference at last investigation were respectively 135 cm, 39 kg, and 51 cm, respectively (Table 1).We Sanger sequencing showed that the disease status segregated with the variants, that is, the variant is present at the homozygous state only in the proband, whereas it is heterozygous in his 2 unaffected sisters and his parents (Figure 1A).
The third individual is affected by mild ID, autism, and attention-deficit/hyperactivity disorder (Table 1).He was born from Iranian consanguineous parents and exome sequencing revealed the presence of a homozygous 4bp-deletion leading to a frameshift in NSUN6 (NM_182543.5:c.1320_1323del;NP_872349.1:p.(Glu441Profs*15)).Sanger sequencing showed the presence of the frameshift at the heterozygous state in the 2 parents and confirmed its homozygous state in the affected son (Figure 1A).This variant is reported only at the heterozygous state in 51, 2, and 62 individuals, respectively, in gnomAD (MAF = 0.0002002), the Iranome database (MAF = 0.00125), and TopMed (MAF = 0.000 4938) and a MAF = 0.00004 in the REGENERON database.Contrary to the frameshift variant of individual 1, the mapping of this frameshift variant in the last exon of NSUN6 suggests that the transcript likely escapes nonsense-mediated decay and is translated into a truncated protein, a hypothesis we could not test because of the unavailability of the patient's cells or RNA.Would our hypothesis be correct, this would lead to a truncated protein lacking the last 29 amino acids in particular Phe458, which, similar to Asp323, is pivotal for the recognition and methylation of the C72 of the tRNA, a known methylation target of NSUN6. 53,55,56Asp323 is essential for the binding of the methyl donor SAM and both NSUN6-Asp323Ala and NSUN6-Phe458Ala lack enzymatic activity. 53Threedimensional protein modeling shows that both Asp323 and Phe458 are positioned in the catalytic site with Asp323 in direct contact with C72 of the tRNA (Figure 1B and C, Supplemental Figure 1).Although the missense change Asp323Asn seems minor, the presence of the amine group in the Asn sidechain does not accommodate the C72 at the same position and the negative charge of Asp, which possibly stabilizes the orientation of the tip of Lys248, would be lost (Figure 1B).

Functional analyses of the identified NSUN6 variants
We engineered the identified missense and the last-exon truncating NSUN6 variant to characterize them further functionally.We tried to express in Escherichia coli an aminoterminally histidine-tagged NSUN6 protein lacking the last 29 residues that mimics the Iranian variant.These encompass the last helix and, notably, the last β-strand (Phe458-Lys465), which is located between β-strands Gly366-Thr372 and Gln395-Gln397 to create an antiparallel β-sheet.Thus, absence of the Phe458-Lys465 C-terminal strand is likely to severely affect the protein fold (Supplemental Figure 2).Consistent with this hypothesis, no soluble protein could be harvested unlike wild type, Asp323Asn, and other missensecontaining NSUN6 proteins produced previously. 53We were similarly unable to purify the same truncated protein tagged at its carboxy-terminal end.Our results confirm that the Cterminal amino acids are important for the overall architecture of NSUN6.
We then assessed the enzymatic activity and the ability to bind the methyl donor SAM of the identified missense variant p.(Asp323Asn).Purified tRNA Cys (GCA) and tRNA Thr (UGU) were incubated with 6 histidine-tagged NSUN6-Asp323Asn and NSUN6-WT and analyzed by methylation assay and RNA mass spectrometry.Whereas m 5 C-modified tRNAs were identified upon incubation with the wild-type enzyme, m 5 C level was undetectable in presence of NSUN6-Asp323Asn (Figure 2A and B).Our results show that the Asp323Asn variant identified in the affected Pakistani individual 2, like the previous assessed engineered mutant NSUN6-Asp323Ala, cannot catalyze m 5 C modification on tRNAs 53 (Figure 2A and B).We then tested the ability of the NSUN6-Asp323Asn variant to bind SAM by ITC and confirmed that Asp323 is essential for SAM binding because NSUN6-Asp323Asn has lost the ability to bind this methyl donor (Figure 2C) in contrast to a high binding affinity of NSUN6-WT to SAM reported previously (K d = 1.7 μM ± 0.1). 53In conclusion, our results confirmed that the function of the assessed variants is impaired.

Drosophila model
We investigated whether the loss of the enzymatic activity of NSUN6 is associated with defects in higher cognitive functions.Toward this aim, we took advantage of the fact that Drosophila melanogaster has a clear one-to-one ortholog of NSUN6, CG11109 (hereafter called Nsun6), which acts in a conserved manner. 49We generated Nsun6 knockout mutants using the CRISPR-Cas9 system and verified the deletion by real-time quantitative reverse transcription polymerase chain reaction.The resulting mutant has a deletion of 585bp in exon 2 that leads to the loss of the Nsun6 (Figure 3A and B).No obvious morphological defects were observed in the Nsun6 mutants compared with the wild-type flies.Nsun6 mutants displayed locomotion defects in the negative geotaxis assay with significantly less flies climbing 10 cm in 10 seconds (Figure 3C).Next, we aimed to test the learning ability of Nsun6 mutants using the courtship conditioning paradigm.Importantly, naive Nsun6 males did not display any defect in their normal courtship behavior (Ci WT : 0.52 vs Ci Nsun6 : 0.62) (Figure 3D).To quantify the learning ability of the flies, we determined the time they spent courting upon entrainment with a nonreceptive female.As expected, both entrained wild-type flies and Nsun6 mutants displayed a significant reduction in the courtship behavior in comparison to the naive males (trained Ci WT : 0.07; trained Ci Nsun6 : 0.32) (Figure 3D), but the learning index of Nsun6 mutant flies was significantly reduced compared with wild-type flies (Figure 3E).Our results indicate that loss of Nsun6 impairs cognitive functions.

Discussion
We present 3 unrelated consanguineous individuals carrying different homozygous deleterious variants in NSUN6.Affected individuals present with a neurodevelopmental disorder characterized by ID, global DD, motor delay, and behavioral anomalies (eg, aggressiveness or attentiondeficit/hyperactivity disorder) (Table 1).Although the truncating variant that maps to the first exon is predicted to lead to the loss of NSUN6, we showed that the missense and the last-exon truncating variant impair the methyltransferase activity of the encoded enzyme by affecting, respectively, its catalytic site or its proper folding.The variable severity of the phenotype could result from differences in the level of NSUN6.For example, we cannot exclude residual enzymatic activity from the last-exon frameshift variant that is associated with impaired folding and the less severe phenotype (individual 3).
NSUN6 best described targets are the C72 of tRNA Cys and tRNA Thr . 47,53,56Impairment of their function would be consistent with the fact that most RNA methyltransferases previously linked to ID target tRNA.However, NSUN6 was recently shown to also target multiple mRNAs at a CTCCA consensus sequence motif. 26,49,50,57Although these modifications do not seem to alter the stability of the targeted transcripts, they are associated with an increase in translation levels. 26,47,50,57Most of the targeted mRNAs encode protein-and RNA-binding proteins of genes implicated in neurodevelopmental disorder (eg, NOVA2 and HUWE1), suggesting that NSUN6 might regulate neurodevelopment through fine-tuning of their expression.Some of the mRNA targets were identified and validated taking advantage of Nsun6 homozygous knockout mouse strains.Although these animals exhibited no apparent phenotypes contrary to our Drosophila strain ablated for Nsun6, they were not assessed for cognitive function, and brain was not within the tissues screened for NSUN6 methylase targets. 26,50This discrepancy warrants further deep phenotyping of both Nsun6 mouse and Drosophila knockouts, as well as a screen for NSUN6 targets in the central nervous system (eg, in genes linked to ID) to better understand the role of this enzyme in brain development.Some NSUN6 paralogs have been associated with syndromic ARID.NSUN2 is linked to a form of ID associated with a wide array of ocular symptoms, brain anomalies, chronic nephritis, and Dubowitz syndrome-like traits, such as facial dysmorphism, microcephaly, and growth retardation, [36][37][38][39][40][41][42][43][44] whereas NSUN3 variants caused oxidative phosphorylation deficiency with lactic acidosis, mitochondrial encephalomyopathy, failure to thrive, cerebral atrophy, and seizures. 28,58Although NSUN5 is hemizygous in patients with Williams-Beuren syndrome, 46 the homozygous ablation of its mouse ortholog affects cognition and brain morphology.0][61][62][63] Systematic assessment of mice homozygously deleted for Nsun5 by the International Mouse Phenotyping Consortium also showed decreased locomotor activity and changes in blood cell composition.These phenotypic differences probably reflect both differences in the expression patterns of NSUN genes 35 and in the targets of the encoded enzymes.For example, NSUN5 methylates ribosomal RNA impairing growth and lifespan, 60,64 whereas NSUN2 and NSUN6 are responsible for the lion's share of m 5 C mRNA modifications.In contrast, NSUN3 modifies the mitochondrial tRNA Met facilitating its binding to the nonuniversal AUA and AUU codons and allowing mitochondrial translation to be efficient. 28,45n conclusion, we add NSUN6 to the list of RNA methyltransferases associated with monogenic ID by identifying 3 families with homozygous deleterious variants in its encoding gene.This new gene-phenotype link will readily help with diagnosis and reproductive options.
c.967G>A; NP_872349.1:p.(Asp323Asn) in affected individual 2 born from consanguineous Pakistani parents.Despite the relatively mild amino acid change (Grantham score = 23), this variant has a Combined Annotation Dependent Depletion score of 32 and is predicted to be deleterious by different prediction tools: PolyPhen2 = probably damaging (1); SIFT = Deleterious (0); MutationTaster = Disease causing (0.99); and GERP score = 5.78.It is present in 7 heterozygous individuals from South Asia in gnomAD but not at the homozygous state (MAF = 0.00002797).It is absent in other gnomADassessed populations, as well as from the TopMed, the Iranome, and the REGENERON Genetics Center databases.

Figure 1
Figure 1 Identification of biallelic variants in NSUN6. A. The pedigree and genotypes of the 3 reported families.On the bottom, the localization of the identified NSUN6 variants at the DNA (exons are depicted in the upper part) and protein level.In the schematic representation of the protein, dark blue indicates the RRM, light blue indicates the catalytic core, and green indicates the PUA domain.B. NSUN6 3D protein modeling.Asp323 is in direct contact with the base 72 of the tRNA in the NSUN6 structure bound to tRNACys (GCA)-G2A:C71U.In the image are shown Asp323, Lys248, and C72.The 3′ end of the tRNA is the end of the orange ribbon on the middle right of the figure.The NSUN6 structure is the gray ribbon; tRNACys (GCA)-G2A:C71U is the orange ribbon.C. The relative positions in the NSUN6 structure bound to tRNACys (GCA)-G2A:C71U of Asp323 and Phe458, respectively, mutated and absent in the affected individuals, are shown.The NSUN6 structure and the tRNACys (GCA)-G2A:C71U are depicted with the same color than in panel B. Asp323 shown in red space-filling model is in direct contact with the base 72 of the tRNA in the NSUN6 structure bound to tRNACys (GCA)-G2A:C71U shown in blue space-filling model.The Phe458 in pink space-filling model that would be missing from a truncated protein lacking the last 29 amino acids and is positioned in the catalytic site is also shown.A zoom-in view from the back of the same model can be seen in Supplemental Figure 1.PUA, pseudouridine synthase and archaeosine transglycosylase; RRM, RNA recognition motif; tRNA, transfer RNA.

Figure 2
Figure 2 Functional of the identified missense variant NM_182543.5:c.967G>A;NP_872349.1:p.(Asp323Asn).A. The methylation assay showing that NSUN6-Asp323Asn could not catalyze m 5 C deposition in transcribed tRNACys (left panel) and tRNAThr (right panel) substrates in vitro.The amount of methyl moiety incorporated into the tRNA is shown on the y-axis over time.NSUN6-WT signal is reported in orange-circle, NSUN6-Asp323Asn in light blue-square, and the negative control in light green-triangle.Error bars represent mean ± SD for 3 independent experiments.B. The m 5 C modification of tRNA substrates after NSUN6 (left panel) and NSUN6-Asp323Asn (right panel) incubation was detected on tRNAThr and tRNACys by mass spectrometry.A m 5 C standard as positive control is reported (top panel).The specific nucleoside to base ion mass of m 5 C is presented in orange, whereas the A one is presented in blue.C. The isothermal titration calorimetry of SAM titrated into hNSUN6-Asp323Asn indicates that the mutated protein has lost the ability to bind SAM.Both the raw data of the heat change measured as differential power (DP) over time are reported (left panel) and the data fitted to a one-site binding model, presented as enthalpy changes (ΔH) over molar ratio of the SAM-binding (right panel) are shown.m 5 C, 5-methylcytosine; SAM, S-adenosyl-L-methionine; tRNA, transfer RNA; WT, wild-type.

FFigure 3
Figure 3 The loss of Nsun6 causes locomotion and learning defects.A. Schematic of the Nsun6 locus and the deletion generated by CRISPR/Cas9.B. qPCR analysis of the Nsun6 transcript level in WT and Nsun6 mutant flies.Rpl15 was used as a normalization control.Bars show median ± s.e.m. of biological duplicates.C. Negative geotaxis assay.Relative number of flies of the indicated genotype that climb over 10 cm in 10 seconds.Bars represent the mean ± s.e.m. of 6 independent measurements with 10 flies each.P values were determined with an unpaired, two-tailed t test.D. Quantification of courtship behavior of naive and trained WT and Nsun6 mutants.The Ci is defined as time spent courting/total time.Bars show median ± s.e.m. of all replicates (WT naive: n = 13, WT trained: n = 14, Nsun6 naive: n = 11, and Nsun6 trained: n = 16).P values were determined with an unpaired, two-tailed t test.E. Quantification of the learning index of WT flies and Nsun6 mutants.The learning index is defined as: (Average of Ci naive -Ci trained)/Average of Ci naive.Bars show median ± s.e.m.P values were determined with an unpaired, two-tailed t test.Ci, courtship index; s.e.m., standard error of the mean; WT, wild-type.