Systematic analysis of snRNA genes reveals frequent RNU2-2 variants in dominant and recessive developmental and epileptic encephalopathies

Abstract

Variants in spliceosomal small nuclear RNA (snRNA) genes RNU4-2 (ReNU syndrome), RNU5B-1, and RNU2-2 have recently been linked to dominant neurodevelopmental disorders (NDDs), revealing a major, previously overlooked role for noncoding snRNAs in human disease. Here, we systematically analysed 200 potentially functional snRNA genes in a French cohort comprising 26,911 individuals with rare disorders and through international collaborations. We identify de novo and biallelic variants in RNU2-2 associated with both dominant and recessive NDDs in 126 individuals from 108 unrelated families. Recessive RNU2-2 NDD is at least twice as frequent as the dominant NDD caused by n.4G>A and n.35A>G, and often arises from a de novo variant in trans with an inherited allele, reflecting the high mutability of snRNA genes. Dominant and recessive RNU2-2-NDDs share overlapping clinical features with frequent epilepsy. Blood transcriptomics and DNA methylation analyses revealed subtle, variant-specific effects on splicing and episignatures. Our findings support a gradient-of-impact model and a continuum between dominant and recessive inheritance, establishing RNU2-2 variants as a frequent cause of NDDs, nearly as prevalent as ReNU syndrome.

Extended Data Fig.1.. Distribution of variants across 200 putatively functional genes in the PFMG cohort.

Supplementary Fig. 1 shows the same distribution in the SeqOIA and Auragen subcohorts separately.

Extended Data Fig. 2.. Variant allele fraction (VAF) distributions of snRNA variants in the aggregated PFMG cohort.

a, VAF distribution of de novo variants in snRNA genes with coverage ≥10. b, VAF distribution of biallelic variants in snRNA genes with coverage ≥10. c, Detailed VAF distribution for snRNA genes harbouring variants in the PFMG cohort. Supplementary Fig. 2 shows the same data in the seqOIA and Auragen subcohorts separately.

Extended Data Fig. 3.. gnomAD allele count (AC) distribution of snRNA variants in the aggregated PFMG cohort.

a, De novo variants with coverage ≥10. b and c, Biallelic variants with coverage ≥10 altogether (b) or stratified according to their inheritance (c). Supplementary Fig. 3 shows the same data in the SeqOIA and Auragen subcohorts separately.

Extended Data Fig.4.. Circos plot depicting preferential associations of biallelic variants.

This scheme suggests that, in affected individuals with compound heterozygote variants, variants preferentially co-occur with one affecting the 5′ domain and the other the 3′ domain, suggesting domain-specific combinatorial effects. Variant colouring is the same as in Fig. 3a.

Extended Data Fig.5.. Examples of biallelic variants identified by genome sequencing in the PFMG cohort.

a, Integrative Genomics Viewer (IGV) screenshots of BAM file alignments showing segregation of the RNU2–2 n.127_118del variant in the homozygous state in three affected siblings, with both parents carrying the variant in heterozygous state. b, Left: IGV screenshot displaying the hemizygous RNU2–2 n.149A>T variant in two affected siblings and their heterozygous parent. Right: Copy-number analysis of genome sequencing data revealing a 642 kb deletion (hg38: chr11:62,781,800–63,424,493) encompassing RNU2–2 and 31 additional genes, present in the two affected siblings and the other parent. c, IGV screenshot showing the presence of the RNU2–2 n.142_153dup variant in homozygous state in the affected child and in heterozygous state in both parents.

Extended Data Fig. 6.. Principal component analysis and classification performance of the RNU2–2 episignature.

a, Hierarchical clustering of the clinical features (n=58, rows) of patients with de novo dominant or biallelic RNU2–2 variants (n=74, columns). Categorical data was converted to 0–1 scale, and values were Z-score scaled for each row. Blue-yellow-red scale depicts Z-scores. Lower values indicate a more favourable phenotype, while higher values represent a more severe phenotype. Missing values are shown in grey. Columns are coloured based on the variant classification: dark blue: de novo, dominant (n.4G>A or n.34A>G); light blue: de novo other (VUS); orange: compound, recessive; purple: homozygous, recessive. b, Principal component analysis of clinical features in RNU2–2 variant carriers. Missing values were imputed as 0. Variant colouring is the same as in Ext Fig. 6a.

Extended Data Fig. 7.. Principal component analysis and classification performance of the RNU2–2 episignature.

In all panels, controls are shown in green, dominant n.4G>A carriers in red, dominant n.35A>G carriers in orange, and recessive biallelic cases in blue. a, PCA of adjusted methylation levels at differentially methylated positions (n=201) up to component 6, after correction for expected baseline methylation level based on age at sampling, sex and estimated blood cell counts (n=92 individuals in total) showing separation of dominant n.4G>A carriers from normal controls on axis 1, and the weaker separation between dominant n.35A>G carriers or recessive biallelic cases from normal controls along axis 6. The percentage of variance explained is provided for each component within the axis title. b, PCA of adjusted methylation levels at differentially methylated positions (n=201), after correction for expected baseline methylation level based on age at sampling, sex and estimated blood cell counts on the restriction to normal controls and carriers of the dominant n.35A>G variant along the first two principal components (n=74 individuals in total). The percentage of variance explained is provided for each axis. c, PCA of adjusted methylation levels at differentially methylated positions (n=201), after correction for expected baseline methylation level based on age at sampling, sex and estimated blood cell counts on the restriction to normal controls and carriers of recessive biallelic variants along the first two principal components (n=78 individuals in total). The percentage of variance explained is provided for each axis. d, Predicted probabilities from a three-block cross-validation using a four-class SVM classifier (control, dominant n.4G>A, dominant n.35A>G, recessive biallelic).

Fig. 1.. Systematic in silico analysis of possible functional snRNAs.

a, Filtering strategy used to retain genes encoding possible functional snRNAs. b, and c, Number and distribution of annotated spliceosomal snRNAs before (b) and after (c) filtering. d, Expression of snRNAs in the human brain from ENCODE small RNA-seq data. Upper panel: all mapped reads (including multi-mapped). Lower panel: uniquely mapped reads. Expression is shown as log₁₀ of the maximum normalized RNA-seq signal. Colours indicate the proportion of the snRNA length covered by mapped reads: 100% (purple), 75–99% (dark blue), 50–74% (green), <50% (yellow).

Fig. 2.. Identification of potential novel snRNA gene–disease associations in the PFMG cohort.

The cohort was divided into solved (n=8,343) and unsolved (n=18,568) cases for discovery analyses, with cases solved by variants in snRNAs with known disease association merged into the unsolved group. We compared the proportion of cases with rare variants (gnomAD allele count < 100) between solved and unsolved groups for rare de novo variants (a) and rare biallelic variants (b). Fisher’s test was used to test statistical enrichment in unsolved versus solved cases for genes in which at least 10 patients had variants (minimum number needed to reach statistical significance in the cohort). P-values are shown above the bars.

Fig. 3.. Overview of RNU2–2 variants identified in this study.

a, Two-dimensional predicted structure of U2–2 snRNA showing structural and functional domains. Arrow heads indicate point variants identified in this study. Variants are coloured according to their inheritance: dark blue: de novo, dominant (n.4G>A or n.34A>G); light blue: de novo other (VUS); orange: compound, recessive; purple: homozygous, recessive. The numbers in black represent the count of patients with each variant, for variants identified in more than one family. Other variant types are shown with dotted (deletions), dashed (duplications) or dotted-dashed (indels) lines. The nucleotide differences between RNU2–2 and RNU2–1 are shown using IUPAC codes. Green numbers refer to the numbering of U2–2 nucleotides. Ψ: pseudouridine, m: 2’-O-methyl residues; m6: N6-methyladenosine; ^2,2,7m3Gppp: 2,2,7-trimethylguanosine cap. Green-shaded regions: functional domains of U2 involved in spliceosomal activity: U2/U6 helix II (nt 1–13); branch-point-interacting stem loop (BSL; nt 25–45); and U2/BS branch helix (nt 32–44). Grey-shaded region: Sm site. b, Locations of RNU2–2 variants on linear (unfolded) U2–2 snRNA, showing clustering of variants in the U2/U6 helix, BSL regions, and Sm site (grey-shaded), as well as preferential associations of compound heterozygous variants. Variants are ordered by inheritance mode and by the position of the first variant on the snRNA. Variant colouring is the same as in Fig. 3a.

Fig. 4.. Variant segregation and facial features associated with RNU2–2 variants.

(Retracted for Med Rxiv submission)

Fig. 5. Possible functional impact of RNU2–2 variants.

Overview of RNU2–2 variants identified in this study. a, Two-dimensional structure of the U2 small nuclear RNA (green). Orange boxes indicate variants from this study, with point mutations and single nucleotide insertions represented on the graph along with their respective changes. ^2,2,7m3Gppp: 2,2,7-trimethylguanosine cap, Ψ: pseudouridine, m: 2’-O-methyl residues. b, Zoom-in box representing the two-dimensional predicted structure of the U2 (green) and U6 (red) small nuclear RNAs (snRNAs) U2/U6 helix II during formation of precatalytic spliceosome. Interactions stabilizing these structures as well as mutations potentially affecting their stability are represented (PDB 5XJC). Shaded region: Lsm/Sm protein binding site of U6 snRNA. c, Zoom-in box representing the two-dimensional predicted structure of U2/U6 helixes Ia and Ib. Nucleotides involved in the structure and associated with pathological variants are represented (PDB 5XJC). d, Close-up of U2/BS branch helix where the two-dimensional predicted structure of the interaction between branch point recognition region of U2 snRNA and the intronic branch site (light brown) is represented. The branch point adenosine is highlighted in purple. ‘YNYURAY’: consensus of branch site in metazoans. The mutations related with this region are depicted on the structure (PDB 5XJC). e, Zoom-in box representing the two-dimensional predicted structure of 17S U2 snRNA and the branch point interacting stem loop (BSL) is highlighted. In blue the adenine at position 30 that interacts with the cytosine (purple) at position 40 linked to two reported variants. In black the nucleotides involved directly on the U2/BS formation. The structure of 17S U2 snRNAs’ BSL structure is represented (PDB 7Q3L). f, Close-up of the U2 snRNA Sm ring (light green) where the nucleotides with the mutations close to or within it are marked in orange (PDB 5XJC).

Fig. 6.. Variant-specific alternative splicing perturbations in RNU2–2 variant carriers.

a, Number of significant alternative splicing events (ΔPSI > 0.05) detected using rMATS and linear regression (p-val<0.01), comparing RNU2–2 individuals (dominant: n.4G>A, n=5; n.35A>G, n=5; recessive: n=4) to 49 controls. Splicing categories are colour-coded: exon skipping (SE, blue), alternative 5′ splice sites (A5SS, orange), alternative 3′ splice sites (A3SS, green), and intron retention (RI, red). b, Venn diagram showing minimal overlap of exon skipping events shared among dominant (n.4G>A, n.35A>G) and recessive (biallelic) variant groups. c, Clustermaps of alternative splicing events for individuals carrying n.4G>A (left, n=5), n.35A>G (middle, n=5), and biallelic variants (right, n=4). d, Principal component analysis based on residuals of PSI values after linear regression showing separation of controls (grey, n=49) from individuals with n.4G>A (n=5, blue), n.35A>G (n=5, green), and biallelic variants (n=4, red). Right panel: PC1 vs. PC2; left panel: PC3 vs. PC4. e, Sashimi plots illustrating an exon-skipping isoform shift in ACLY observed in RNU2–2 variant carriers compared to controls.

Fig. 7.. Methylation profiles across differentially methylated probes in RNU2–2 variant carriers and normal controls.

In all panels, controls are shown in green, dominant n.4G>A carriers in red, dominant n.35A>G carriers in orange, and recessive biallelic cases in blue. A, Heatmap of adjusted methylation levels at differentially methylated CpGs. Selected probes are represented in rows while patients and normal controls are represented in columns after hierarchical clustering of their methylation profiles. Probes are grouped according to their association pattern: probes common to all variants, probes specific to n.35A>G, probes specific to biallelic recessive variants, and probes specific to n.4G>A. Samples are annotated by type (Control, Dominant n.4G>A, Dominant n.35A>G, Recessive biallelic) but sample clustering was unsupervised and blind to this classification. b, Raincloud plots showing the average methylation level per sample within each probe association pattern. Normal controls and variant carriers display significantly different distributions, with p-values from two-sided Wilcoxon rank-sum tests indicated above the plots.