De novo variants in the non-coding spliceosomal snRNA gene RNU4-2 are a frequent cause of syndromic neurodevelopmental disorders

A highly recurrent insertion explains 0.52% of undiagnosed NDD in Genomics England

We identified a highly recurrent single base insertion (GRCh38:chr12:120,291,839:T:TA; n.64_65insT) in RNU4-2 in GEL¹. This variant was initially identified as arising de novo in 38 probands recruited for genome sequencing with their unaffected parents¹³. Extending the search to include probands without data for both parents in the full GEL cohort, we identified an additional eight individuals with the n.64_65insT variant; in all eight, the detectable inheritance is consistent with the variant having arisen de novo (i.e. where a single parent sample was available the variant was not detected in it). All of the 46 individuals with the variant have undiagnosed NDD (categorised as global developmental delay, intellectual disability, and/or autism spectrum disorder), corresponding to 0.52% of 8,841 probands with currently undiagnosed NDD in GEL. The n.64_65insT variant is not found in any of 3,408 NDD probands with an existing genetic diagnosis, 21,817 probands with non-NDD phenotypes, or in 33,122 unaffected individuals. Individuals with the variant are significantly enriched for global developmental delay (n=37; OR=3.56; Fisher’s P=2.75×10^-4), delayed gross motor development (n=26; OR=2.55; P=1.64×10^-3), microcephaly (n=26; OR=6.62; P=7.87×10^-10), delayed fine motor development (n=24; OR=2.61; P=1.69×10^-3), hypotonia (n=18; OR=3.60; P=7.09×10^-5), short stature (n=15; OR=3.54; P=2.17×10^-4), drooling (n=7; OR=19.2; P=2.83×10^-7), and absent speech (n=6; OR=6.23; P=7.45×10^-4) compared to all other probands with NDD in GEL (n=12,203; diagnosed and undiagnosed) (Figure 1A; Supplementary Table 1).

• Download figure
• Open in new tab

Figure 1: Characterisation of individuals with the n.64_65insT variant in GEL.

(A) The proportion of individuals with human phenotype ontology (HPO) terms corresponding to phenotypes observed in ≥ 5 individuals with the n.64_65insT variant compared to all other individuals with NDD. Terms that are significantly enriched in individuals with the n.64_65insT variant are marked with a *. Multiple terms relating to global developmental delay, intellectual disability, hypotonia, seizure, microcephaly, autism, and short stature have been collapsed into single phenotypes. Of note, this figure relates only to HPO terms entered for each individual into GEL, which may be incomplete. A more detailed phenotypic characterisation of individuals with variants in RNU4-2 is provided below. (B-D) Quality control metrics for the variant calls in all 46 individuals with the variant: (B) genotype quality scores, (C) allele balance, and (D) coverage.

The n.64_65insT variant is not found in 76,215 genome-sequenced individuals in gnomADv4.0¹⁴, or in 245,400 individuals in the All of Us dataset¹⁵. It is seen in a single individual in the UK Biobank¹⁶ (allele frequency=1.02×10^-6) with a variant allele balance consistent with a true variant (23 reference and 18 [44%] alternate reads). This individual has an ICD-10 code for ‘personal history of disease of the nervous system and sense organs’ but no further phenotype data to assess a potential NDD diagnosis (Supplementary Table 2).

Given the high recurrence rate of this insertion, we wanted to rule out that it is a sequencing or mapping error, despite the overwhelming evidence of phenotype enrichment. Notably, the variant is a single A insertion after a run of four Ts, ruling out the most common cause of sequencing error for indels, polymerase slippage in homopolymer repeats. The variant calls were all high quality based on both analysis of quality metrics (Figure 1B-D) and manual inspection on IGV (Supplementary Figure 1). Finally, the genomic region surrounding the insertion and RNU4-2 maps uniquely to a single region of the genome with short-read sequencing in GRCh38 and T2T CHM13v2.0/hs1 (Supplementary Figure 2).

The n.64_65insT variant is within a highly constrained region with multiple NDD-causing variants

The recurrent n.64_65insT variant resides within the central region of RNU4-2, towards the 5’ end of an 18 bp region which is depleted of variants in population datasets compared with the rest of the gene (26% of all possible SNVs observed in UK Biobank compared to a median of 78% across the rest of the gene; Figure 2A; Supplementary Figure 3). Based on the population variant data, we defined a critical, highly constrained region as chr12:120,291,825-120,291,842.

• Download figure
• Open in new tab

Figure 2: A highly structured 18 bp region of RNU4-2 that is critical for BRR2 helicase activity is enriched for variants in NDD and depleted in population cohorts.

(A) Allele frequency of variants in 7,519 undiagnosed NDD probands GEL (teal) and the UK Biobank cohort (grey) across RNU4-2. The 18 bp critical region is marked by a horizontal bar at the top of the plot. (B) Schematic of U4 (teal) binding to U6 snRNA (grey). The 18 bp critical region is underlined. (C) The structure of U4 and U6 snRNAs resolved by cryoEM²¹. Created using RCSB Protein Data Bank²² (structure 6QW6). In both (B) and (C) single base insertions identified in individuals with NDD are shown by black arrows and positions of SNVs by red nucleotides.

We searched for variants across this region in GEL, and also in additional cohorts containing undiagnosed individuals with NDD (see methods). In total, we identified 119 individuals with variants across this region (Table 1), the vast majority of which have the initial n.64_65insT variant (n=92; 77.3%). For 92 of the 119 individuals, sequencing data for both parents was available to confirm the variants had arisen de novo. Five of the 11 additional variants are also single base insertions, including n.77_78insT (GRCh38:chr12:120,291,826:T:TA), which is seen in six individuals, two of whom are affected siblings. The enrichment of single base insertion variants in this region in individuals with NDD is striking: 54/8,841 (0.61%) GEL undiagnosed NDD probands (55/10,388 individuals) have single base insertions compared to 2/490,132 individuals in the UK Biobank (OR=1,531; 95%CI:404,>16,384; Fisher’s P=3.3×10^-92).

Aside from insertions, there is also a modest enrichment of SNVs in GEL NDD probands across the critical region (undiagnosed NDD: 6/8,841; UK Biobank 35/490,132; OR=9.51; 95%CI:3.27-22.8; Fisher’s P=8.16×10^-5). We identified 16 individuals across cohorts with SNVs in this region (Table 1; 11 confirmed de novo), all with phenotypes consistent with individuals with insertion variants. The identified SNVs cluster with the two regions harbouring insertion variants at the extreme ends of the 18 bp critical region (Figure 2B). Conversely, SNVs in the central portion (particularly at nucleotides 71-74) are observed in both non-NDD individuals in GEL (n=2) and population controls, although all at low frequencies (Supplementary Table 3). Across the remainder of RNU4-2 there is no significant enrichment of variants in undiagnosed NDD probands when compared to non-NDD probands (194/7,519 undiagnosed NDD; 521/19,428 non-NDD in GEL aggregated variant dataset²⁰; OR=0.96; 95%CI:0.81-1.14; Fisher’s P=0.67).

In total, we identify variants in this 18 bp region in 119 individuals with NDD. This includes 60/8,841, or 0.68%, of all genetically undiagnosed NDD probands in GEL (0.49% of all NDD probands). In contrast, variants in this region are observed in 39/490,132 (0.008%) individuals in the UK Biobank (OR=85.8; 95%CI:56.4-131.6; Fisher’s P=1.84×10^-78).

U4 snRNA binds to U6 snRNA through extensive complementary base-pairing in the U4/U6.U5 tri-snRNP complex of the major spliceosome. Unwinding of U4 and U6 is essential to generate the catalytically active spliceosome². The 18 bp critical region in RNU4-2 maps to a single-stranded region of U4 between the stem I region of complementary base-pairing to U6 and the 3’ stem-loop structures (nucleotides 62 to 79; Figure 2B). This region is known to be loaded into the active site of the SNRNP200-encoded BRR2 helicase, which mediates unwinding of the U4/U6 duplex². The highly recurrent n.64_65insT variant is within a previously described ‘quasi pseudoknot’, or T-loop, structure²¹ (Figure 2C). The region spanning nucleotides 76 to 78, where the recurrent n.77_78insT variant resides, is involved in base-pairing with U6 in stem III²³ (Figure 2C). Both of these regions are thought to stabilise the U4/U6 interaction and accurately position the U6 ACAGAGA sequence to receive the 5’splice site during spliceosome activation. Insertion of a single base into either of these structures may destabilise the U4/U6 interaction and/or alter the positioning of the U6 ACAGAGA sequence and potentially disrupt the correct loading of the 5’ splice site into the fully assembled spliceosome. Nearby regions that are predicted to have important roles, such as the U4/U6 stem I binding region, are not enriched for variants in NDD probands.

Variants in this crucial region cause a severe syndromic NDD phenotype

To characterise the phenotypic spectrum associated with variants in RNU4-2, we collected detailed phenotypic information for a subset of 36 individuals (33 with n.64_65insT, one with n.64_65insG, and two with SNVs; Table 2; Supplementary Table 4). Using these data, we find the RNU4-2 syndromic NDD to be characterised by moderate to severe global developmental delay (two children with SNVs with moderate delay) and intellectual disability in all individuals. The majority (82%) achieved ambulation but at a delayed age (average 3.6 years, range 18 months to 7.5 years) with many noted to have a wide-based or ataxic gait. Only one individual (with an SNV) had fluent speech, some had a few words, and most were non-verbal. All but one were reported to have dysmorphic facial features. These facial features varied but consisted of a myopathic face with deep set eyes (some widely spaced and some narrowly spaced), epicanthus, wide nasal bridge, anteverted nares or underdeveloped ala nasi coll, large cupped ears (some posteriorly rotated), full cheeks, a distinctive mouth with full lips with downturned corners, high arched palate, and a large or protruding tongue.

Associated growth and neurodevelopmental phenotypes present in ≥75% of individuals include short stature, microcephaly (mostly congenital), seizures (spanning infantile spasms, focal seizures and generalised tonic-clonic seizures, febrile seizures, and status epilepticus with variable onset from the first year of life, but most between 3-10 years of age), and hypotonia. Brain MRI showed a spectrum of abnormalities in the majority of individuals, most frequently reduced white matter volume, non-specific abnormalities of the white matter, hypoplasia of the corpus callosum, ventriculomegaly, and delayed myelination. Involvement of multiple organ systems was reported for all individuals, often including visual (optic nerve hypoplasia, cortical blindness, strabismus, nystagmus), gastrointestinal (constipation, reflux, feeding issues with need for a gastrostomy tube), and bone/skeletal abnormalities (osteopenia, recurrent fractures, scoliosis, kyphosis, hip dysplasia), and in a lesser number of individuals, hearing, endocrine (hypothyroidism, growth hormone deficiency), limb, sleep, genitourinary, dental, cardiac, and cutaneous concerns (Table 2; Supplementary Table 4).

RNA-sequencing in blood does not show a global disruption to splicing

Given the importance of U4 snRNA in the spliceosome and previous observations of global disruption to splicing observed in other spliceosomopathies²⁴, we analysed RNA sequencing data from blood samples for five individuals from GEL. Three of these individuals have the highly recurrent n.64_65insT variant, another has the other recurrent insertion, n.77_78insT, and the final patient has an SNV (n.78A>C). We did not see any significant difference in the number of gene expression outliers using OUTRIDER²⁵, or in the number of retained introns, or all outlier events using FRASER2²⁶ in the five individuals with RNU4-2 variants compared with 5,409 controls (Supplementary Table 5). At present, RNA from additional tissues (e.g. brain samples) of affected individuals is not available. It is possible that the observed RNU4-2 variants disrupt more subtle aspects of alternative splicing in a tissue-specific manner, as has been observed for other snRNA variants²⁷.

RNU4-2 is highly expressed across tissues and in the brain across development

Humans have multiple genes that encode the U4 snRNA, although only two of these, RNU4-2 and RNU4-1, are highly expressed in the human brain (Supplementary Table 6). RNU4-2 and RNU4-1 are contiguous on chr12, both 141 bp long, and highly homologous, differing by four nucleotides (97.2% homology). RNU4-1 has a similar depletion of variants in population cohorts in the centre of the RNA, however, we do not observe an enrichment of variants in GEL in this central region. Further, we do not observe an enrichment of variants in GEL in this central region (Supplementary Figure 4). There is a variant equivalent to our highly recurrent variant in RNU4-1 that is observed in six individuals in the UK Biobank dataset. There are no consistent phenotypes recorded in these six individuals (Supplementary Table 2).

To investigate the reason for variants in RNU4-2, but not RNU4-1, causing NDD, we analysed the expression of both RNU4-1 and RNU4-2 in the brain. First, we analysed the expression patterns of both genes across multiple developmental stages using bulk RNA-seq data from 176 human prefrontal cortex samples in BrainVar²⁸. The expression of RNU4-1 and RNU4-2 is tightly correlated (Supplementary Figure 5), however, RNU4-2 is consistently expressed at a significantly higher level than RNU4-1 (Figure 3A). Secondly, we assessed chromatin accessibility in the chromosome 12 locus containing both RNU4-1 and RNU4-2 using ATAC-seq data from two human prenatal prefrontal cortex samples. These data show a dramatic chromatin accessibility signal around RNU4-2 and a much lower signal surrounding RNU4-1, again consistent with much higher expression of RNU4-2 in the brain (Figure 3B). Overall, these data support the role of RNU4-2 as the major U4 transcript in the brain.

• Download figure
• Open in new tab

Figure 3: RNU4-2 is more highly expressed than RNU4-1 in the prefrontal cortex.

(A) Levels of RNU4-1 (grey) and RNU4-2 (teal) expression at different developmental stages from BrainVar²⁸. (B) ATAC-seq data from human prenatal prefrontal cortex (18 and 19 gestational weeks (GW)) with substantially higher peaks of chromatin accessibility around RNU4-2 (teal) than RNU4-1 (grey).

Multiple factors likely explain the high recurrence of the n.64_65insT variant

The n.64_65insT variant is highly recurrent. It is observed in 46/12,249 NDD probands in GEL (0.38%; or 0.52% of undiagnosed NDD probands). In contrast, the most recurrent protein-coding variant in a dataset of 31,058 individuals with developmental disorders²⁹ is observed in 36 individuals (0.12%; GRCh38:chr11:66211206:C:T; PACS1:p.Arg203Trp). The exact reasons for this high recurrence are unclear, however, we hypothesise three contributing factors. First, a high local mutation rate, which may be driven by the open chromatin state and very high levels of transcription (Figure 3). In UK Biobank, a median of 76% of all possible SNVs in RNU4-2 are observed (calculated across 18 bp sliding windows). This is compared with 13% on average in 1,000 random intergenic sequences of the same length (141 bp; P<0.001, Monte-Carlo Fisher-Pitman test; Supplementary Figure 6). Despite the high number of variants in RNU4-2 in UK Biobank, there are no individuals with homozygous variants and all observed variants are very rare (maximum allele frequency = 0.025%), consistent with high levels of selection acting on variants across RNU4-2.

Secondly, a high overall mutational burden does not explain the high recurrence of this specific single base insertion. Local formation of secondary structure and base stacking is a known driver of biased small insertion mutations³⁰. The high propensity of this region to form secondary structure when single-stranded may drive creation of this specific insertion. Finally, it may be that germline selection is acting to increase the frequency of this specific variant, as has been shown for other highly recurrent sites³¹. While we see no association with paternal age (mean 33.1 in probands with RNU4-2 variants and 33.4 across other NDD probands; Supplementary Figure 7), fully testing this hypothesis will require deep sequencing of testes or sperm samples.

No other spliceosomal snRNA genes are enriched for de novo variants in NDD

Given the newly identified importance of RNU4-2 in NDD, we sought to determine whether other snRNA genes with no known association to NDD could also harbour novel diagnoses. We investigated 28 snRNA genes that are expressed in the brain, using multiple approaches (Supplementary Table 7). First, we tested for an overall enrichment of de novo variants in undiagnosed NDD probands compared to non-NDD probands across each snRNA with at least two identified de novo variants in probands with undiagnosed NDD (n=13) using the high-confidence de novo callset in GEL. Of the 12 genes other than RNU4-2, none showed a significant enrichment of de novo variants in undiagnosed NDD probands (all Fisher’s P>0.15).

Secondly, hypothesising that the burden of pathogenic variants in other snRNAs may be restricted to specific critical regions, as we see for RNU4-2, we used an 18 bp sliding window to identify snRNA regions that are depleted of variation in the UK Biobank compared to the overall variant burden across each gene. Notably, the regions with the highest depletion in RNU4ATAC correspond to two hotspots of pathogenic variants in ClinVar (chr2:121530923-121530946, chr2:121530984-121531007), however, the strength of the depletion in these regions is lower than in RNU4-2 (minimum normalised proportion of observed −0.11 and −0.2 versus −0.5 for the depleted region in RNU4-2), consistent with lower selection acting on variants in RNU4ATAC that cause recessive disorders. We identified 14 regions in 13 unique snRNAs with a deviation from the median number of SNVs across the full gene of at least 20% (Figure 4; Supplementary Table 8). We repeated our de novo variant enrichment test in regions with at least two de novo variants in undiagnosed NDD probands (n=3). Only the conserved region in RNU4-2 was significant (Fisher’s P=1.34×10^-9; undiagnosed NDD probands n=33, non-NDD probands n=0; all other tests Fisher’s P>0.25).

• Download figure
• Open in new tab

Figure 4: Multiple snRNA genes have regions that are depleted of variation in the population.

The proportion of observed SNVs in 490,640 genome sequenced individuals in the UK Biobank, in sliding windows of 18 bp across each snRNA gene, normalised to the median value for each gene.

Finally, we looked for recurrent de novo variants in undiagnosed GEL NDD probands that were absent from diagnosed NDD probands, non-NDD probands, and population controls. There are three de novo variants with an allele count ≥3 in the GEL undiagnosed NDD cohort, two in RNU1-2 (chr1:16,895,992:C:T and chr1:16,896,002:A:G), and one in RNVU1-7 (chr1:148,038,767:G:A). However, all three variants are observed at comparable frequencies in non-NDD probands and are also found at relatively high frequencies in population controls (all variants’ AF>0.5% in gnomAD 4.0).

Categorising participants in Genomics England

We defined four groups of individuals in GEL v18. Individuals with NDD (n=13,812) were defined as those with human phenotype ontology (HPO)³⁵ and/or International Classification of Diseases 10th Revision (ICD-10) codes³⁶ for global developmental delay (HP:0001263, HP:0012736, HP:0011344, HP:0011343, HP:0011342; ICD-10: R62, F80, F81, F82, F83, F88, F89), intellectual disability (HPO: HP:0001249, HP:0002187, HP:0010864, HP:0002342, HP:0001256, HP:0006887, HP:0006889; ICD-10: F70, F71, F72, F73, F78, F79), and/or autism (HPO: HP:0000717, HP:0000729, HP:0000753; ICD-10: F84), or who were recruited to GEL with a normalised specific disease of intellectual disability. NDD individuals were classified as diagnosed (n=3,424) if they were marked as solved or partially solved in the gmc_exit_questionnaire table or had an entry in the submitted_diagnostic_discovery table in GEL Labkey. The remaining 10,388 NDD individuals formed our undiagnosed NDD cohort. Of these, 8,841 are probands. We also identified 21,817 probands without NDD phenotypes (i.e. without the HPO and ICD10 codes detailed above) and 33,122 individuals reported to be unaffected. Our defined cohorts exclude anyone who has subsequently removed consent.

For the majority of our analyses, we used two previously defined datasets within GEL. First, a high-confidence set of de novo variants from 13,949 trios¹³. As of 13 March 2024, this dataset includes 12,554 probands with consent: 5,426 probands with undiagnosed NDD, 2,352 with diagnosed NDD, and 4,776 non-NDD probands. De novo variants were filtered to those that pass the stringent_filter. Second, an aggregated variant call set (aggV2)²⁰ which contains 29,850 probands: 7,519 undiagnosed NDD, 2,903 diagnosed NDD, and 19,428 non-NDD.

Identifying variants in population datasets

We used data from gnomAD v4.0 (76,215 genome sequenced individuals)¹⁴, All of Us¹⁵ (accessed via the publicly available data browser https://databrowser.researchallofus.org/; 245,400 genomes as of 28 March 2023) and the UK Biobank (490,640 genome sequenced individuals)¹⁶.

Expanded NDD cohort and clinical data collection

Clinical data were collected from research participants after obtaining written informed consent from the parents or legal guardians, with the study approved by the local regulatory authority. Samples were collected largely through personal communications (NW, AODL, DGM) as variants in this gene have not been prioritised in analysis. On entry into Matchmaker Exchange using the seqr node, one match was made (CD). NW reviewed the National Health Service Genome Medicine Service (NHS GMS; V3) dataset. Samples from NHS GMS were manually checked to remove duplicates with GEL. AODL reviewed the Broad Center for Mendelian Genomics and the GREGoR consortium datasets. DGM contacted additional local collaborators. Clinical data were collected and summarised for features seen across the cohort.

We additionally searched 7,149 trios with autism spectrum disorder and 4,180 sibling control trios from three cohorts: Simons Simplex Collection (SSC; 2,383 cases; 1,938 controls)¹⁸, SPARK (3,144 cases; 2,190 controls)³⁷, and MSSNG (1,622 cases; 52 controls)¹⁷.

Generating 1,000 random intergenic sequences

Using the bedtools subtractBed function³⁸ we retrieved regions on chromosome 12 that do not overlap with RefSeq transcripts aligned by NCBI. We further removed regions within 10 kbp of an annotated transcript and restricted the remaining regions to those at least 141 bp in length (n=611). We further removed regions overlapping the centromere. We then generated a set of 1,000 random sequences from each intergenic region and then randomly selected 1,000 non-overlapping regions from these.

Identifying human snRNA genes

We extracted genes with snRNA biotypes from Ensembl genome annotation v111. We filtered out known pseudogenes (i.e. with gene names marked with “P” or identified through manual curation). For each remaining gene, we used BrainVar²⁸ RNA-seq expression data to calculate the mean CPM value across the gene. We selected only genes with mean CPM value across all BrainVar samples >5, resulting in a dataset of 28 snRNA genes.

Assessing variant depletion

Given the high mutability of RNU4-2 and other snRNA genes, coupled with strong selection pressures on variants, we did not think that conventional mutational models would be well calibrated to assess variant depletion. Instead, we devised a sliding window-based strategy to identify regions within snRNA genes that are relatively depleted of SNVs. We split genes into 18 bp sliding windows (chosen as it is the size of the region defined in RNU4-2) and tallied the number of SNVs observed in UK Biobank 500k genome sequencing data within that window, divided by the total number of possible SNVs (i.e. 18x3). The proportion of possible SNVs observed in each window was normalised to the median across all sliding windows in that gene (i.e. the per-gene median proportion observed was subtracted from each value). Depleted regions were defined as those spanning windows with a deviation from the per-gene median of at least 20%, i.e. normalised observed proportion of possible SNVs < −0.2. The same calculation was performed on 1,000 randomly selected 141 bp intergenic regions on chr12 (see above). A one-way approximative (Monte Carlo) Fisher-Pitman test was conducted to show the median observed proportion of possible SNVs was significantly higher for RNU4-1 and RNU4-2 compared to the distribution in the 1,000 random regions.

RNA-sequencing of individuals with RNU4-2 variants

Blood was collected from a subset of 100,000 Genomes Project probands in PaxGene tubes to preserve RNA at the time of recruitment. RNA was extracted, depleted of globin and ribosomal RNAs, and subjected to sequencing by Illumina using 100 bp paired-end reads, with a mean of 102M mapped reads per individual. Alignment was performed using Illumina’s DRAGEN pipeline. FRASER2²⁶ and OUTRIDER²⁵ were used to detect abnormal splicing events and expression differences with samples run in batches of 500, both run via the DROP pipeline³⁹.

Analysing RNU4-2 and RNU4-1 expression

We used the BrainVar²⁸ dataset to assess patterns of whole-gene expression of RNU4-2 and RNU4-1 in the human cortex across prenatal and postnatal development. This dataset includes bulk-tissue RNA-seq data from 176 de-identified postmortem samples of the dorsolateral prefrontal cortex (DLPFC, n=167 older than 10 post-conception weeks) or frontal cerebral wall (n=9 younger than 10 post-conception weeks), ranging from 6 post-conception weeks to 20 years of age. The 100 bp paired-read RNA-seq data from BrainVar were aligned to the GRCh38.p12 human genome using STAR aligner⁴⁰, and gene-level read counts for GENCODE v31 human gene definitions were calculated with DEXSeq⁴¹ and normalised to counts per million (CPM)⁴².

Prenatal prefrontal cortex ATAC-seq data

Methods of generating ATAC-seq have been described previously⁴³, which is the source of the data shown here. Briefly, fresh prenatal (18 and 19 gestational weeks) brain samples were dissected within 2 hours of elective termination to extract the entire telencephalic wall, from the ventricular zone to the meninges. Intact nuclei were isolated by manually douncing the tissue on ice using a loose pestle douncer then lysed on ice for 10 minutes by adding a solution with 0.1% NP-40. Nuclei were spun down by centrifugation then resuspended and exposed to Tagmentation Enzyme for 30 minutes at 37C. The ATAC-seq library was generated using Illumina barcode oligos, amplified by high-fidelity PCR, and sequenced on the Illumina HiSeq 2500 using paired-end sequencing. Reads were aligned to GRCh38 using the ENCODE ATAC-seq pipeline with default parameters⁴⁴. A UCSC Browser track of per nucleotide ATAC-seq counts was used to assess the region around RNU4-2 and RNU4-1.

Burden testing and statistical analysis

The enrichment of de novo variants across each of 28 snRNA genes and 14 constrained sub-regions was assessed in undiagnosed NDD probands compared to non-NDD probands using the high-confidence de novo callset. Odds ratios and associated P-values were calculated using Fisher’s exact test in R. A P-value threshold of 0.0031 was used to assess statistical significance as a Bonferroni correction accounting for 16 tests.