A non-syndromic orofacial cleft risk locus links tRNA splicing defects to neural crest cell pathologies [RNA-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. Molecularly, 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, and impaired 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, and the concomitant accumulation of tRNA fragments, affect neural crest and craniofacial development and positioned MYCN, DDX1, and tRNA processing defects as risk factors in the pathogenesis of orofacial clefts. Overall design: Bulk-RNA-sequencing experiments were performed in triplicates for DDX1 KO cNCCs cultured with 0.4 ng/mL DOX (DDX1 ON) and without DOX for 4 days (DDX1 OFF). cDNA libraries were prepared using the RNA sample preparation kit (TruSeq v2; Ilumina) as previously described (Respuela et al., 2016) and sequenced with a 2x 100-bp strand-specific protocol on a HiSeq 2500 sequencer (Illumina). Basic read quality control was performed using FastQC (Babraham Bioinformatics) and read statistics were obtained with SAMtools. Reads (61-73 million per sample) were mapped to the human reference assembly (GRCh38.p13) using STAR Aligner (Dobin et al., 2013). Aligned reads were sorted with featureCounts (Liao et al., 2014), specifying for paired end reads. Differential expression analysis was performed using DESeq2 (Love et al., 2014), excluding low-count transcripts (<5).

口面裂（Orofacial clefts）是全球范围内最常见的先天性颅面畸形（congenital craniofacial malformations）。此类出生缺陷的病因呈多因素性，涉及遗传与环境因素。但在多数病例中，其潜在致病机制仍未明确，阻碍了对该疾病分子发病机制的解析。
本研究整合了全基因组关联数据（genome-wide association data）、病例与对照队列的靶向重测序（targeted re-sequencing）、细胞类型特异性表观基因组谱分析（cell type-specific epigenomic profiling）以及基因组结构分析（genome architecture analyses），从功能与分子层面解析了与非综合征型口面裂风险升高相关的一个基因组位点。

我们发现，与口面裂相关的常见及罕见风险变异，与一个保守增强子（conserved enhancer）e2p24.2存在交集——该保守增强子在颅神经嵴细胞（cranial neural crest cells）中被激活，而颅神经嵴细胞是负责塑造颅面部结构的胚胎细胞类型。我们将e2p24.2的远程基因组相互作用定位至一个包含MYCN与DDX1基因的拓扑关联结构域（topologically associated domain），并证实MYCN与DDX1均为鸡胚颅面部发育所必需。

分子水平上，我们发现e2p24.2在颅神经嵴细胞中调控MYCN的表达，但对DDX1并无此调控作用。而DDX1是MYC家族转录因子（transcription factors）的靶基因，同时也是tRNA剪接复合物（tRNA splicing complex）的组成成分。颅神经嵴细胞中DDX1的缺失会导致未剪接tRNA片段（unspliced tRNA fragments）的积累，并同时损伤全局蛋白质合成（global protein synthesis）与颅神经嵴细胞的迁移能力。我们进一步证实，诱导tRNA片段的产生足以破坏颅面部发育。

综上，本研究揭示了一条分子机制：tRNA剪接受损伴随tRNA片段积累会影响神经嵴细胞与颅面部发育，并将MYCN、DDX1及tRNA加工缺陷确定为口面裂发病的风险因素。

总体实验设计：
本研究针对两种培养条件下的DDX1敲除颅神经嵴细胞（cranial neural crest cells，简称cNCCs）开展了三次生物学重复的批量RNA测序（Bulk-RNA-sequencing）：一组添加0.4 ng/mL多西环素（Doxycycline，DOX）以维持DDX1表达（DDX1 ON组），另一组不添加DOX以关闭DDX1表达（DDX1 OFF组），处理时长均为4天。
参照既往研究（Respuela等，2016）的方法，使用RNA样本制备试剂盒（TruSeq v2；Illumina）构建cDNA文库（cDNA libraries），并采用链特异性2×100 bp测序策略在HiSeq 2500测序仪（Illumina）上完成测序。
使用FastQC（巴布拉汉生物信息学研究所，Babraham Bioinformatics）完成基础测序读段质量控制，通过SAMtools获取读段统计信息。使用STAR序列比对工具（STAR Aligner）将每个样本的6100万至7300万条读段比对至人类参考基因组组装版本GRCh38.p13。使用featureCounts工具对双端读段的比对结果进行排序。使用DESeq2工具进行差异表达分析，过滤掉计数低于5的低表达转录本。