Massively parallel reporter perturbation assays uncover temporal regulatory architecture during neural differentiation

Abstract

Gene regulatory elements play a key role in orchestrating gene expression during cellular differentiation, but what determines their function over time remains largely unknown. Here, we perform perturbation-based massively parallel reporter assays at seven early time points of neural differentiation to systematically characterize how regulatory elements and motifs within them guide cellular differentiation. By perturbing over 2,000 putative DNA binding motifs in active regulatory regions, we delineate four categories of functional elements, and observe that activity direction is mostly determined by the sequence itself, while the magnitude of effect depends on the cellular environment. We also find that fine-tuning transcription rates is often achieved by a combined activity of adjacent activating and repressing elements. Our work provides a blueprint for the sequence components needed to induce different transcriptional patterns in general and specifically during neural differentiation.

Fig. 1. Experimental design.

a Computational framework to select regions and perturbation sites. b Heatmap of motif instances in the assayed regions (left); distribution of the number of motifs perturbed in each region (top right); distribution of the number of regions harboring each motif (bottom right). c Library design; selected regions were included in their wild-type (WT) form, selected motifs were perturbed (by altering the sequence in the predicted motif site) using three perturbation methods individually (PERT single) as well as in combination with other perturbations in selected cases (PERT double). Random sites were perturbed (RAND) and the entire WT sequence was scrambled as negative controls (SCRAM) for each WT sequence. d The designed sequences were synthesized and cloned into the lentiMPRA vector and associated with 15-bp barcodes. ARE antirepressor element, BC barcode. Reporter, EGFP enhanced green fluorescent protein, LTR long terminal repeat, mP minimal promoter, WPRE Woodchuck Hepatitis Virus Posttranscriptional Regulatory Element. e lentiMPRA libraries were infected into hESCs and following 3 days, we induced neural differentiation via dual-SMAD inhibition and obtained DNA and RNA at seven time points (0, 3, 6, 12, 24, 48, and 72 h). f Association between barcodes and designed sequences, and the number of barcodes observed in DNA and RNA sequencing was determined using MPRAflow. Differential analysis between WT and PERT activity to determine motif regulatory effect over time was assessed using MPRAnalyze. Source data are provided as a Source Data file.

Fig. 2. Preprocessing, consistency, and categorization of FRSs.

a Illustration of the four filters applied to perturbed sequences to remove inactive and non-functional sites, both at each timepoint and across timepoints (error bars represent mean ± 1 SD; “Methods”). b Number of sequences that passed all three filters for each perturbation method. c Definition of main and sub-categories of motif binding effects based on their effect on transcription. d Distribution of categories across FRSs that pass the four filters and are under the same main category (activating or dampening) in perturbation method 3 and at least one of the perturbation methods 1 or 2. The distribution is shown across FRSs (top) and across the unique motifs and regions composing the FRSs in this study (bottom).

Fig. 3. Examples of the four sub-categories of FRSs.

In each figure, the WT sequence is indicated in black, centered at the mean activity and including error bars of ±1 SD across the three replicates, and SCRAM in blue including error bars of mean ± 1 SD across all scrambled sequences. Each motif is plotted in a different color, including error bars of mean ± 1 SD across the three replicates, and its perturbation effect in the regions is indicated in the text box. a NEOROD2 has a silencer effect. b NANOG, TP73, SOX1 have contributing, inhibiting, contributing effects concordantly. c BX088580.2, TP73 both have essential effects. All genomic coordinates are hg19. Source data are provided as a Source Data file.

Fig. 4. Temporal motif effects.

The a activating or b repressing motifs in at least one time point. red—activator motifs, blue—repressor motifs. Color scale indicates the average of the perturbation signal (LogFC(WT/PERT)) across all significant instances of motifs for a specific TF (row normalized). Data are organized using hierarchical clustering and early, mid, and late clusters are indicated. Genome browser snapshots of assaysed sequences near predicted sequence motifs that are associated with c OTX2 and d BARHL1 TFs, showing the motifs that were perturbed and their effect on activity across time points. e HOXD9 repressor motif example. Line plots similar to Fig. 3, mean activity $\pm$ 1 SD across the three replicates. All coordinates are hg19. Source data are provided as a Source Data file.

Fig. 5. Double perturbation scheme.

a Experimental design for perturbing two single motifs separately and then a double perturbation of both simultaneously, and the requirements for being included in downstream analysis. b Volcano plot for the model testing for log-additivity of the individual effects. c–g Examples of double perturbation results demonstrating different patterns of cooperation: log-additive effects consistent with a billboard model (c); log-additive effects of one dampening and one activating element (d); fully dependent cooperation consistent with the enhanceosome model (e); a billboard-enhanceosome hybrid model with one required element and one with a dampening effect (f); a redundancy example, perturbing either motif has negligible effect, but perturbing both has a substantial effect (g). All coordinates are hg19. Line plots similar to Fig. 3, mean activity ± 1 SD across the three replicates. Source data are provided as a Source Data file.