A non-syndromic orofacial cleft risk locus links tRNA splicing defects to neural crest cell pathologies [Ribo-seq]

Orofacial clefts are the most common form of congenital craniofacial malformations worldwide. The etiology of these birth defects is multifactorial, involving genetic and environmental factors. In most cases, however, the underlying causes remain unexplained, precluding molecular understanding of disease mechanisms. Here, we integrated genome-wide association data, targeted re-sequencing of case and control cohorts, cell type-specific epigenomic profiling, and genome architecture analyses, to functionally and molecularly dissect a genomic locus associated with an increased risk of non-syndromic orofacial cleft. We found that common and rare risk variants associated with orofacial cleft intersect with a conserved enhancer (e2p24.2) that becomes activated in cranial neural crest cells—the embryonic cell type responsible for sculpting the craniofacial complex. We mapped e2p24.2 long-range interactions to a topologically associated domain harboring MYCN and DDX1 and demonstrated that both MYCN and DDX1 are required for craniofacial development in chicken embryos. We found that e2p24.2 regulates the expression of MYCN, but not DDX1, in cranial neural crest cells. In turn, DDX1 is a target of the MYC family of transcription factors and a component of the tRNA splicing complex. The loss of DDX1 in cranial neural crest cells resulted in the accumulation of unspliced tRNA fragments, depletion of the mature pool of intron-containing tRNAs, and ribosome stalling at codons decoded by these tRNAs. These effects were accompanied by defects in both global protein synthesis and cranial neural crest cell migration. We further showed that the induction of tRNA fragments is sufficient to disrupt craniofacial development. Together, these results uncovered a molecular mechanism in which impaired tRNA splicing affects neural crest and craniofacial development and positioned MYCN, DDX1, and tRNA processing defects as risk factors in the pathogenesis of orofacial clefts. Overall design: Ribosome profiling was performed as described (McGlincy and Ingolia, 2017) with minor modifications. cNCCs derived from three separate differentiations of DDX1-/-::TetON-HA-DDX1 hiPSCs were grown for 5 days with (DDX1-ON) or without (DDX1-OFF) DOX to ~80% confluency on 10 cm dishes, rapidly washed twice with ice-cold PBS supplemented with 100 µg/mL cycloheximide and lysed on-plate in 500 µL ice-cold Polysome Lysis Buffer (20mM HEPES pH 7.3, 100mM KCl, 5mM MgCl, 1% TritonX-100, 2 mM DTT, 100 µg/mL cycloheximide). Lysate was passed through a 26 gauge needle 10 times and clarified by centrifugation for 5 minutes at 1,300 x g, 4°C, and the supernatant was further clarified by centrifugation for 5 minutes at 20,000 x g, 4°C. Samples were normalized to the same concentration, as measured by NanoDrop Spectrophotometer, and treated with 350 U Ambion RNAseI (Invitrogen) per O.D. of absorbance at 260 nM per mL for 16 hours at 4°C. Samples were loaded on cold 10-50% sucrose gradients generated by layering 6 mL of 10% sucrose gradient buffer (20 mM HEPES pH 7.3, 100 mM KCl, 5 mM MgCl, 1 mM DTT, 100 µg/mL cycloheximide, 10% sucrose w/v) on 6 mL of 50% sucrose gradient buffer (20 mM HEPES pH 7.3, 100 mM KCl, 5 mM MgCl, 1 mM DTT, 100 µg/mL cycloheximide, 50% sucrose w/v) and mixed using a BioComp Gradient Station (BioComp, New Brunswick, Canada). Sample-loaded gradients were spun at 210,000 x g for 2.5 hours at 4°C using a SW-41 swinging-bucket rotor (Beckman Coulter) and Optima XPN Ultracentrifuge (Beckman Coulter). These gradients were fractionated on a BioComp Gradient Station and the monosome peaks were collected. RNA was extracted from the collected monosome peaks using Trizol LS (Invitrogen) and purified using a Direct-zol RNA Miniprep kit (Zymo Research). Fragments between 17-34 nt long were isolated and libraries were prepared from them as described in McGlincy and Ingolia, 2017, except with a depletion of rRNA immediately before reverse transcription using a mixture of rRNA-complementary biotinylated oligos (http://bartellab.wi.mit.edu/protocols.html). Reads were processed for downstream analysis as follows: Basic read quality control was performed using FastQC (Babraham Bioinformatics). Illumina libraries were demultiplexed, adapter sequences were trimmed using fastp (Chen, 2023). Libraries were demultiplexed by using UMI-tools to assign correctable barcodes to reads, and split into samples using SeqKit2 (Shen et al., n.d.; Smith et al., 2017). Reads aligning to rRNA were identified with BowTie and removed (Langmead et al., 2009). The remaining reads were aligned to the genome using human reference assembly (GRCh38.p13) with STAR (Dobin et al., 2013). STAR was run using the following relevant parameters: --outFilterType BySJout --outFilterMismatchNmax 2 --outFilterMultimapNmax 1 --outFilterMatchNmin 16 --alignEndsType EndToEnd. Aligned reads were then deduplicated using the 7 nt unique molecular identifier with UMI-tools (Smith et al., 2017). Read lengths and their phasing were quantified on maximal spanning windows generated using plastid (Dunn and Weissman, 2016). Differential analysis was performed by counting the ribosome profiling reads aligned to the CDS of each gene with featureCounts (Liao et al., 2014). Differential expression analysis was performed on counts using DESeq2 (Love et al., 2014), excluding genes with fewer than one read in all samples. The DESeq2 results can be found in Supplementary Table 1. Principal component analysis was performed on regularized log-transformed count data using the rlog() and prcomp() functions in R (R Core Team, 2021). Codon analysis was performed as described previously with minor modifications (Darnell et al., 2018). Aligned reads were trimmed by 12 nt from both ends, and read coverage was calculated across each 300 nt window centered at each in-frame codon of all annotated CDSs in the transcriptome, excluding the first and last 150 nt of each CDS and those with less than 0.33 reads per nt. Read coverage at each position in each codon window was normalized to the mean coverage of their CDS of origin, and these ribosome densities at every position were averaged across all instances for each codon. These positional means were additionally normalized to the total summed ribosome density for each sample before subtracting densities of DDX1OFF samples from DDX1ON samples for each replicate to generate ? ribosome densities.