Computational Inference Tuned with Wet‐bench Results Yields a Model of the Transcriptional Regulatory Network Governing Zebrafish Periderm Differentiation

Abstract

Early embryonic epidermis and oral epithelium in mammals stratify to generate outer protective layer, called the periderm, which prevents water loss and interepithelial adhesions among limbs and oral structures. Orofacial cleft is a structural birth defect which results from the mutations in transcription factors regulating oral periderm differentiation. In zebrafish, enveloping layer (EVL), also called periderm, is a simple squamous epithelial monolayer which arises from blastomeres. EVL differentiation shares many regulatory transcription factors (TFs) with mammalian periderm differentiation including Irf6 and Grhl3. Therefore, we use zebrafish EVL as a tractable model for the difficult‐to‐access mammalian oral periderm. We previously carried out ATAC‐seq in isolated periderm and flow‐through cells and found that the Grhl3, Klf17, Tfap2a, Cebpb, and Gata3 binding sites (TFBS) were enriched in periderm‐specific ATAC‐seq peaks, implicating these TFs in the transcriptional regulatory network (TRN) driving the periderm differentiation. In this study, we applied a computational method to infer the structure of the zebrafish periderm TRN, and, for use in tuning the parameters of the algorithm, we generated a reduced‐representation network model based on biologically‐validated edges. We applied a network‐modeling algorithm called “modified least absolute shrinkage and selection operator with stability approach to regularization selection” (mLASSO‐StARS) to publicly‐available zebrafish single‐cell RNA‐seq datasets and to “prior” evidence of TF‐gene interactions based on the presence of TFBS in enhancers associated with their proximal genes. To generate a reduced‐representation network comprised of biologically‐validated edges we performed RNA‐seq on zebrafish wild‐type and irf6 maternal mutant embryos at 6 hours post fertilization (hpf). Next, in an addback paradigm, we injected into irf6 maternal mutant embryos at the one cell stage mRNA encoding transcription factors Irf6, Grhl3, Klf17, Tfap2a, Cebpb, and Gata3 (separately), and at 6 hpf harvested RNA and profiled expression of 94 genes using nanoString probes. Partially overlapping sets of genes were rescued by each TF. To identify the direct targets of Grhl3, IRF6 and Tfap2a, we performed ChIP‐seq and CUT&RUN‐seq. Results from the addback paradigm and CUT&RUN assays were compiled to generate a gold standard network for periderm differentiation. Finally, we used gold standard network to test the consequence of varying input parameters on the performance of computationally‐derived TRN. Satisfyingly, adding the prior information to the co‐expression network improved the TRN model significantly. The resulting TRN provide a better understanding of periderm differentiation in zebrafish, will identify key hub genes in the network, which in turn may help in prioritization of candidates identified in genomic analyses of orofacial cleft patients.