Parallel genetics of regulatory sequences in vivo

Understanding how regulatory sequences control gene expression is fundamental to explain how phenotypes arise in health and disease. Traditional reporter assays inform about function of individual regulatory elements, typically in isolation. However, regulatory elements must ultimately be understood by perturbing them within their genomic environment and developmental- or tissue-specific contexts. This is technically challenging; therefore, few regulatory elements have been characterized in vivo. Here, we used inducible Cas9 and multiplexed guide RNAs to create hundreds of mutations in enhancers/promoters and 3′ UTRs of 16 genes in C. elegans. To quantify the consequences of mutations on expression, we developed a targeted RNA sequencing strategy across hundreds of mutant animals. We were also able to systematically and quantitatively assign fitness cost to mutations. Finally, we identified and characterized sequence elements that strongly regulate phenotypic traits. Our approach enables highly parallelized, functional analysis of regulatory sequences in vivo.


Introduction
Understanding gene regulation is fundamental for understanding development and tissue function in health and disease. Animal genomes contain diverse regulatory sequences, which are organized in contiguous stretches of genomic DNA, ranging from a few to hundreds or thousands of bases. Promoters, enhancers, and silencers act mainly on transcription while transcribed sequences (including 5′ and 3′ UTRs) regulate splicing, export, localization, degradation, and translation of mRNAs. Many regulatory sequences encode multiple regulatory functions which can cooperate, compensate and antagonize each other [1][2][3] . Understanding this logic requires combinatory perturbations. Moreover, a single binding site, due to fuzzy recognition motifs, may tolerate certain mutations [4][5][6] . The interaction of eff ectors with regulatory elements can additionally be modulated by sequence structure 2,5,7 , co-factors 2,3 , chemical modifi cations and the temporal order of binding 1,3 . Sequence activity is therefore dependent on native sequence context, cell type, tissue, development and the environment. Accordingly, phenotypic consequences of mutations in regulatory regions are diffi cult to predict. To understand biological functions, an approach to target regulatory sequences in vivo with many diff erent mutations is required.
Although parallel interrogation of regulatory sequences has been developed in cell lines and yeast [8][9][10][11][12] , only a few in vivo approaches have been achieved in animal models. These use integration of reporters 13,14 or injection of RNA libraries 15,16 and therefore do not evaluate endogenous phenotypes, or are restricted to one stage of the animal life cycle. Classical genome editing by injection, now widely accessible due to CRISPR-Cas 17,18 , has enabled endogenous functional tests, but is still work-intensive and limited in scalability 19,20 .
Here we use inducible expression of Cas9 and multiplexed single guide RNAs in C. elegans populations to generate hundreds of targeted mutations in parallel. We targeted diff erent regulatory regions across 16 genes and analyzed >12,000 Cas9-induced mutations to fi rst describe characteristics of dsDNA break repair in the C. elegans germline and the genotype diversity that this introduces at the targeted loci. We then applied our mutagenesis approach to generate hundreds of deletions along the well-studied lin-41 3′ UTR, which is targeted by the miRNA let-7 [21][22][23][24] . We developed an RNA-sequencing-based strategy to quantify the impact of each mutation on lin-41 RNA level. Using DNA sequencing we followed the relative abundance of these diff erent mutations over several generations to infer their phenotype. This analysis revealed a previously undescribed compensatory function between the two let-7 miRNA binding sites. Finally, we couple the targeted mutagenesis of regulatory sequences to selection of phenotypic traits. We isolate 57 reduction-of-function alleles in 3 genes that show strong changes in phenotype, mediated by mutations in enhancer, TATA-box, and 3′ UTRs which act by both loss and gain of gene regulation.

Inducible Cas9 to generate hundreds of mutants in parallel.
To introduce many diff erent, targeted mutations in vivo, we developed an approach in C. elegans using inducible expression of Cas9 and several multiplexed "single guide RNAs" (sgRNAs). This required only few initial injections to create transgenic animals, allowed maintenance without mutagenesis, and enabled time-controlled creation of mutated populations in parallel, with sizes only limited by culturing approaches (up to ~10 7 in our case). Mutant populations could then be used for various purposes. For example, they could be selected by phenotype or reporter activity, or directly analyzed by targeted sequencing to measure the impact of mutations on RNA regulation or fi tness (Fig. 1a).
As an initial test, we generated transgenic lines with plasmids for heat shock-driven Cas9 expression and oneor multiple sgRNAs targeting a ubiquitously expressed single-copy GFP reporter. After a transient heat shock, we could observe GFP-negative animals in culture, indicating activity of Cas9. We performed a two-hour heat shock induction of Cas9 in the parents (P0) and collected progeny (F1) in a time course experiment. The highest fractions of mutants were obtained 14 -16 hours after heat shock, with ~50% (sg1) and ~20% (sg2) of eggs producing GFP negative animals (Supplementary Fig. 1a). We obtained similar results when targeting the dpy-10 gene and counting  Fig.1b).
Characteristic CRISPR-Cas9-induced mutations from 91 GFP negative animals consisted of deletions, insertions or a combination of both and originated from sgRNA cut sites (Supplementary Fig. 1c, d). When we used three sgRNAs within the same transgenic line targeting adjacent positions, deletions appeared around one cut site or spanned between two cut sites (Supplementary Fig. 1d). This indicated that pools of sgRNAs could lead to more diverse genotypes and cover more nucleotides. Most deletions induced by a single sgRNA were between 3 -10 bp long and we observed insertion lengths between 1 -30 bp (Supplementary Fig.  1e).
Homozygous animals will be produced in the F2 by heterozygous self-fertilizing F1. Additionally, since Cas9 induced in the P0 could still be active after fertilization, F1 animals could be mosaic with a wild type germline and mutant somatic cells (Supplementary Fig. 1f). We therefore wanted to assess how many germline mutations were generated. For this we analyzed the inheritance of GFP negative animals from F1 to F2 generations using an automated fl ow system and found that ~80% of mutations were indeed germline mutations (Supplementary Fig. 1gi). For the rest of our work we used such non-mosaic F2, generated by F1 germline mutations.
To analyze large mutated C. elegans populations in bulk, we established a targeted sequencing protocol based on long 0.5-3 kb PCR amplicons. It allowed us to sequence the complete targeted locus, to capture large deletions and to multiplex samples for sequencing (Supplementary Fig.  2a). This resulted in ~200,000 -800,000-fold read coverage at the targeted regions. We also created a software pipeline ("crispr-DART" for "CRISPR-Cas Downstream Analysis and Reporting Tool") to handle targeted sequencing data of such amplicons and analyze the contained indels. The pipeline works with various targeted sequencing technologies and extracts and quantifi es indels. The output contains html reports of coverage, mutation profi les, sgRNA effi ciencies and optional comparisons between two samples. Processed genomics fi les from the output can then be used for more in-depth custom analyses with additionally supplied R scripts (Supplementary Fig. 2b-d) ("Code Availability" in Supplementary Methods).
To test our approach in larger scale, we induced Cas9 in 50,000 P0 animals by heat shock, and amplicon sequenced the mutated locus from bulk samples of 400,000 F2 progeny. Deletions per genomic base-pair peaked sharply around sgRNA cut sites. Pools of multiplexed sgRNA plasmids resulted in deletions spanning two or several sgRNAs ("multi cut") in addition to smaller deletions surrounding single sgRNAs ("single cut") ( Fig. 1b). Insertions occurred within a few nucleotides to cut sites and were less frequent than deletions (~1/2 -1/10 th ) (Fig. 1c).

Fig. 1 | Cas9 induction to produce multiplexed C. elegans mutants. a,
Outline of our approach which uses heat-shock Cas9 induction to create large "diversifi ed" populations containing indel mutations at the targeted region. Mutated populations can be used for various downstream assays: selection according to morphological traits and reporter activity, bulk RNA sequencing to measure eff ects of individual 3′ UTR mutations, or DNA sampling over several generations to infer fi tness of diff erent genotypes. b, Example of a targeted locus (upstream of snb-1). The percentage of DNA sequencing reads containing deletions with respect to the total read coverage is plotted at the corresponding genomic position. Bulk worm samples were sequenced, thus 2% deletions per genomic nucleotide refers to approximately 2% of worms with a deletion at the respective nucleotide. SgRNA cut sites are indicated by orange triangles. Individual deletion events are shown below in red. c, Same analysis as in b) but for insertion events.
Finally, we assessed the diversity of the generated mutations. We started by counting the number of unique deletions per base pair. We fi rst studied deletions created by single-cut events for each sgRNA and found that highly active sgRNAs could generate up to 150 unique deletions (rows in Fig. 2e). Most deletions covered a region 10 -12 bp surrounding the cut sites. On average, every sgRNA could generate 15 unique deletions per bp at the center of the cut site and up to an additional 5 unique deletions 5 bp away from the cut site (black line profi le in Fig. 2e). We then studied multi-cut events. Naturally, these deletions spanned much larger regions, as defi ned by the spacing between two or more sgRNAs. Here we found up to 200 unique deletions per base pair and on average 20 unique deletions per sgRNA covering a region >500 bp surrounding each cut site (Fig.  2f). We then considered each unique indel a genotype (>0.001% frequency and >5 reads). On average, one sgRNA created 50 deletion genotypes and 10 insertion genotypes. However, some sgRNAs created up to 400 genotypes (Fig.  2g). Since we used several sgRNAs per transgenic line, we observed on average ~200 insertion and deletion genotypes per sample and in effi cient lines up to 1700 deletion and 1200 insertion genotypes (Supplementary Fig. 5a). More effi cient sgRNAs resulted in a higher number of unique deletions (Supplementary Fig. 5b). Transgenic lines expressing more sgRNAs showed more unique deletion genotypes, possibly because of an increased chance of containing effi cient sgRNAs and the combined activity of multiple sgRNAs (Supplementary Fig. 5c). Together these data show that inducible expression of Cas9 with multiplexed sgRNAs can induce a variety of mutations to study hundreds of regulatory variants in parallel. This includes small deletions to target individual elements at nucleotide resolution, large deletions to interrogate combinatory interactions and insertions to change spacing between sites or duplicate existing sequences.
Quantifying mutant gene expression and fi tness. A major challenge to understand gene regulation is the interaction of diff erent elements. Especially in 3′ UTRs, which can act on all levels of gene expression, this can be diffi cult. To simultaneously measure mRNA levels for all generated 3′ UTR deletions within large C. elegans populations, we developed a targeted RNA sequencing strategy. As a proof of principle, we tested it on a microRNA-regulated mRNA. The lin-41 mRNA is regulated by let-7 microRNAs which bind two complementary sites in the 1.1 kb long 3′ UTR (LCS1 and LCS2, 22 and 20 nucleotides long, separated by a 27 nt spacer) 21,22,24,32,33 (Supplementary Fig. 6a).
Although studies with reporter plasmids showed that each binding site could not function on its own 22 , other studies concluded that each site could recapitulate wild-type regulation when present in three copies 23 . Prior to our work, compensation had not been tested for each site, in the native sequence context, or at natural expression levels. Therefore, we targeted the lin-41 3′ UTR with a pool of 8 sgRNAs and two sgRNA pairs close to the LCSs (Fig. 3a). Lin-41 downregulation occurs with let-7 expression in the developmental stages L3-L4 21,[32][33][34] . To measure let-7-dependent regulation, we collected RNA from bulk worms at L1 and L4 stages. We extracted L4 stage RNA after complete lin-41 mRNA downregulation by let-7 35 and before the occurrence of the lethal vulva bursting phenotype 24 (for an overview of samples see  Supplementary Table 1). a, Effi ciency measured for each sgRNA (n=127) per experiment. This is calculated by counting the fraction of reads containing deletions +/-5 bp around the sgRNA cut site. b, Proportions of reads with diff erent types of mutations detected in each experiment (n=60 experiments). Complex: reads with more than one insertion or deletion, or additional substitutions which suggest a combination of multiple events. c,d, Length distribution of mutations supported by >5 reads and >0.001% frequency, found in all experiments. Deletions were classifi ed as multi cut deletions (n=2915) when a deletion overlapped with more than one sgRNA cut site +/-5bp or otherwise were classifi ed as single cut deletions (n=3169) (for a scheme describing this see Supplementary Fig. 2  . We then sequenced lin-41specifi c cDNA with long reads covering the complete 3′ UTR (Supplementary Fig. 6c). Each read contained full information on any deletion present in the RNA molecule, while the number of reads supporting each deletion could be used to estimate RNA expression level. To determine let-7dependent eff ects, we then analyzed how diff erent deletions aff ect RNA abundance at L4-relative to L1 stage. We observed an average 4-fold up-regulation of lin-41 at larvae stage L4 only when both let-7 seed sites were mutated by deletions (Fig. 3b). No signifi cant eff ect was observed from deletions aff ecting either LCS1 or LCS2 seeds alone. The 4-fold regulatory eff ect is consistent with the known magnitude of down-regulation in the natural context 22,32,33 or the up-regulation when disrupting both let-7 interactions 24,36,37 (2 -4-fold). A very slight and non-signifi cant up-regulation was observed for LCS1. Interestingly, LCS1 has a longer seed pairing (8 nt) than LCS2 (6 nt and a G-U pair) and was detected with more reads in an in vivo miRNA proximity ligation approach 38 . In addition 3xLCS1 acted stronger than 3xLCS2 in a reporter assay 23 . As an independent approach, we used UMAP on the k-mer content of long cDNA reads to obtain clusters of genotypes 39 . These data also suggest that RNA molecules transcribed from clusters with large deletions overlapping both LCS seeds were detected with more reads in L4 stage compared to L1 stage animals (Supplementary Fig. 6d-f).
To assign fi tness to individual mutations in a controlled environment, we established measurements on genotype abundance over several generations. For this, we sampled genomic DNA of consecutive generations (Fig. 3c). Disrupting let-7 regulation of lin-41 is known to result in lethal developmental defects 21,24,32,37,40 . Deletions in the lin-41 3′ UTR which overlapped with both LCS seeds quickly disappeared from the population and were mostly absent already after one generation (Fig. 3d). Consistent with our observed absence of an eff ect on RNA expression, deletions aff ecting either one of the two LCS alone were not signifi cantly depleted. Our results indicate that the two let-7 complementary sites can largely compensate each other's loss under laboratory conditions. We conclude that parallel mutagenesis coupled with targeted RNA or DNA sequencing can be used to analyze function and interactions of regulatory elements in vivo directly from large populations in bulk.
Gene regulatory sequences which aff ect phenotypic traits. Next, we w anted to isolate regulatory sequence variants which aff ected animal phenotype. This could be useful to discover functional elements, provide starting points to study mechanisms, and to obtain animals with desired phenotypic traits. Such an approach would also capture any functional sequences regardless of the level, time or tissue of regulation, which can be hard to predict. We proceeded to select worms based on strong organismal phenotypes aff ecting animal movement and body shape (Unc, Rol, Dpy) and to identify the causative mutations. We targeted a predicted enhancer 41 , three promoters and all 3′ UTRs of 8 genes with known loss-of-function phenotypes and screened 35,000 -45,000 animals for each of these regions (Supplementary Table 2). To determine which mutations were initially present in the screened population, we performed targeted sequencing on siblings (Supplementary Fig. 8a, b). Initially, we isolated several mutants with large deletions (>500 bp) that disrupted the coding frame or the polyadenylation signal (AATAAA) (Supplementary Fig. 8c, d). Similar large-scale, on-target deletions have previously been described in cell lines 42,43 and mice 44 . We also found large insertions (<250 bp) which originated from within +/-1 kb of the targeted region, or from loci on other chromosomes (Supplementary Fig. 8c,  d). We found such large deletions or insertions in 5 out of 8 screened genes, demonstrating that for these genes our screen was sensitive enough to detect animals with aff ected phenotypes (Supplementary Table 2). From the screen we isolated 57 reduction-of-function alleles in 3 genes (egl-30, sqt-2, sqt-3) and none from the other 5 genes (dpy-2, dpy-10, rol-6, unc-26, unc-54) (Supplementary Table 2). Deletions, insertions and complex mutations were represented equally among isolates (Supplementary Fig. 8e). The observed phenotypic traits showed complete penetrance and we scored their expressivity which diff ered between mutations. We found that several mutations in the 3′ UTR of egl-30 lead to Uncoordinated (Unc) phenotypes. In 7/11 mutants, a region circa 100 bp downstream of the STOP codon was aff ected. The smallest deletion was 6 bp (Fig. 4a) (Supplementary  Fig. 8g). We also found mutations with Roller (Rol) phenotypes overlapping a putative sqt-2 enhancer predicted from chromatin accessibility profi ling 41 (Supplementary  Fig. 8f). We also targeted sqt-3, a gene associated with three distinct morphological traits 45,46 (Dpy, Rol and Lon). In total we isolated 39 alleles. 13 mutations upstream of sqt-3 likely aff ected transcriptional initiation, with 11/13 overlapping the predicted TATA-box (Fig. 4b).
In line with the Rol phenotype, which indicates reduction-of-function, pre-mRNA and mRNA levels were both reduced to ~50% (Supplementary Fig. 9d). The remaining sqt-3 alleles were 3′ UTR mutations. Almost all (25/26) were insertions or complex mutations originating at sg2 (Fig. 4c). The only deletion overlapped with the polyadenylation signal (AAUAAA). We knew from sequencing siblings that sg2 was very effi cient (~25%) and that in our samples various deletions covering the 3′ UTR were present. We therefore isolated 13 distinct non-Rol mutants using direct PCR screening ( Supplementary Fig. 9a, b). Despite containing deletions or insertions originating at the effi cient sg2, these animals showed the wild type non-Rol trait (Fig. 4d). We did follow-up experiments with one of the 25 insertion alleles, sqt-3(ins) and determined that mRNA levels were    Fig. 9c, d). Since deletions in this region were well tolerated (non-Rol), we concluded that the isolated Rol mutations likely resulted from a gain of repressive sequence which led to the observed reduction of mRNA. Interestingly the polyA mutant sqt-3(polyA), for which mRNA levels were equally reduced to 50%, showed a much weaker Rol phenotype, with only slight bending of the head (Fig. 4c,  Supplementary Fig 9d, e). This suggests that additional mechanisms besides mRNA down-regulation might reduce protein output in sqt-3(ins).
To defi ne the repressive sequence elements, we targeted the inserted sequence with several sgRNAs and screened for non-Rol "revertants" in which the wild type trait was restored. 12/13 revertants contained deletions overlapping with the insertion, with the smallest being 5 bp (Fig. 4d). A restored wild type trait likely resulted from restored expression levels and indeed mRNA levels in two independent revertants were restored to normal (Supplementary Fig. 10b). Overall, predicted RNA secondary structures did not change, suggesting other factors cause the Rol phenotype of sqt-3(ins) (Supplementary Fig. 10c). Finally, we found, using sequence transplantations into hypodermally expressed dpy-10 and neuronally expressed unc-22, that repression was dependent on 3′ UTR context and observable also in neuronal tissue (Supplementary Fig. 10d). This implies that we isolated a general repressive/destabilizing RNA element.To discover interacting 3′ UTR elements we had included sgRNAs for the remaining 3′ UTR. This revealed a compensatory deletion upstream of the insertion, which was able to revert the Rol phenotype. We isolated two more additional alleles after using sgRNAs specifi c for this region (Fig. 4d) (Supplementary Fig. 10a). Surprisingly, mRNA levels were not restored (Supplementary Fig. 10b). This points to an alternative mechanism of restored protein function, for example on translational level, or aff ecting mRNA at a diff erent developmental time point.
Overall, these results demonstrate that parallel genetics and selection by phenotype can be used to obtain specifi c phenotypic traits, to fi nd functional sequences, and to discover unexpected regulatory interactions in vivo.

Discussion
In this study we develop a general approach for parallel genetics of regulatory sequences in vivo, using inducible expression of a CRISPR-nuclease and multiplexed sgRNAs. Large "diversifi ed" populations can then be used for comprehensive analysis using direct deep sequencing or for selection by phenotypic traits. This allows directly linking regulatory genotypes with phenotypes. We demonstrate this in the model organism C. elegans but believe it is similarly applicable in other animals or plants which allow transgenesis and inducible expression of genome editors.
As we show, sgRNA effi ciencies around 1.5 % are suffi cient to analyze eff ects of mutations on gene regulation and phenotype when coupled with deep sequencing and manual-or automated selection of animals from large populations. However, higher effi ciency would be desired for improved comprehensive testing. This could be achieved with more advanced expression systems 47 and optimization of sgRNA promoters for high germline expression. Alternative induction systems (e.g. Auxin, FLP/FRT, Cre/lox, Gal/UAS) could enable continuous and germline-specifi c Cas9 expression to increase effi ciency and allow directed evolution experiments. Our targeted sequencing protocol can capture long deletions, uses the same amplicon for the whole locus, and allows sample multiplexing. Unique molecule counting methods 48 for long reads 49,50 could be incorporated to reduce PCR biases. Established protocols 28 are available for shorter (100-300 bp) target regions. We assumed that each animal in bulk samples contributed equally to the extracted genomic DNA. In the future, barcoding to determine genotypes of individual animals could be developed with plate-based or split-pool 51,52 methods. Our method only works at nucleotide resolution close to the sgRNA cut sites. To allow denser tiling of regions with mutations, CRISPR-nucleases with alternative or dispensable PAM requirements could be used 17,53,54 . Although indels are applicable to regulatory regions and even coding sequences 55-57 , point mutagenesis would enable fi ne mapping of regulatory nucleotides and amino acids. This exciting possibility could be opened up by implementing hyperactive base 58-61 -or prime editors 62-64 . Alternatives to CRISPR-Cas could be developed based on inducible expression of cassette exchange 65,66 or targeted transposons 67,68 . Targeting several independent loci in one step might be applied to candidate regulatory elements (e.g. miRNA targets, enhancers), screening candidate genes (e.g. from networks, pathways 69 , other assays 70 ), and for synthetic co-evolution of several loci 71 .
Indel data from high throughput genome editing in human cell lines led to insights into dsDNA break repair outcomes [27][28][29][30][72][73][74][75] . Whether these data align with outcomes in vivo or in the germline has not been well studied. We found longer indels and less 1nt templated insertions in our data. This can be explained with higher activity of microhomology-mediated end joining (MMEJ), which has been shown to be the main dsDNA break repair pathway in C. elegans 31 . Mutations typical for MMEJ have been implicated in diseases 76 and microhomologies allow predicting the outcomes of CRISPR-Cas9 editing [27][28][29][72][73][74][75] . Deep mutation profi les in C. elegans could be used to study MMEJ in vivo.
In the evolution and determination of phenotypes gain or loss of regulatory sequences is important 1,77,78 and can be modeled as a gradual process by single nucleotide changes 6,77-81 . Although indels occur less often than single nucleotide mutations naturally, their impact can be more severe. Deletions have already been found to destroy or generate regulatory elements during evolution 77,78 . Insertions templated from regions around dsDNA breaks during DNA repair could duplicate functional sequence elements and might be an underestimated force in the evolution of regulatory sequences or disease mutations.
In line with previous observations we found that the two let-7 complementary sites repress lin-41 ~4-fold 22,24,32,33,36,37 . Furthermore, we uncovered a previously undescribed compensatory interaction. Previous studies concluded either that one LCS could not recapitulate wild-type repression at all 22 or tested only multiple copies of each site 23,82 . These studies were done with reporter plasmids and therefore not in the native sequence context (which aff ects target site accessibility 23 ), or at natural expression levels (which aff ects miRNA-target ratios 37 ). We found a slight but nonsignifi cant eff ect on RNA regulation when disrupting LCS1. Together with previous studies 23,38 this indicates that LCS1 could act more strongly than LCS2, and that both sites cannot completely compensate each other's loss. Since our method might lack the sensitivity to detect more subtle eff ects from deletion of each LCS alone, future studies would be needed to test this thoroughly.
Our screen for sequences that aff ect phenotypic traits doubles the phenotypic regulatory alleles registered in the last forty years at Wormbase. Due to the mutagenesis effi ciency and because we screened for strong phenotypes, we did not saturate and likely missed many mutations. We expect that higher effi ciencies would likely lead to mutations aff ecting expression and phenotype also for the genes for which we did not isolate any mutations. To determine comprehensively which mutations are tolerated by a locus, even higher effi ciencies would be needed. Some mutations were likely aff ecting known regulatory sequences (sqt-3 TATA-box), while for others (sqt-2 enh, sqt-3 and egl-30 3′ UTR) we did not determine a mechanism. However, such mutants can be a starting point for future studies on regulatory mechanisms using biochemical, computational and genetic approaches. Additionally, the approach described here can be applied to isolate alleles with altered function (reduction-of-function, gain-of-function) and to obtain animals with desired phenotypic traits. Furthermore we discovered compensatory deletions in the sqt-3(ins) 3′ UTR that restored the wild type phenotype. Although not restoring mRNA levels, and likely acting by an alternative mechanism, we uncovered these by screening for the animal phenotype. Our results highlight the possibility to uncover unexpected regulatory interactions using in vivo interrogation and readout of phenotype. Such interactions should be detectable, even if they act across diff erent developmental time points, tissues or cell types.
Altogether we believe the approach presented here will help in understanding how gene regulatory logic aff ects phenotype in vivo.

Maintenance of animals.
C. elegans were maintained on NGM plates with Escherichia coli OP50 as originally described 83  sgRNA design. sgRNAs were designed manually using C. elegans genome version ce11 and the plasmid editor Ape (A plasmid Editor, M.W. Davis, https://jorgensen.biology.utah.edu/wayned/ape/) and evaluated using the E-CRISP web application (http://www.e-crisp.org/E-CRISP) 90 . For regulatory regions of interest, we aimed at a regular spacing between target sites, dense coverage and as little as possible predicted off -targets with less than three mismatches. A detailed list of sgRNA sequences, together with their characteristics, effi ciency prediction scores and predicted off -targets, can be found in a supplementary table.
Generation of transgenic C. elegans. Extra-chromosomal array transgenes were generated by standard procedure using micro-injection into the gonad 91 .
A detailed list of injection mixes and their composition can be found in a supplementary table. The injection mix usually contained plasmids for heatshock inducible Cas9 (pMB67 87 or pJJF152) at 50 ng/μL, 1-10 sgRNAs (using the backbone pMB70 87 or pJJR50 88 ) at 10-50 ng/μL, a visual coinjection marker expressing mCherry in the pharynx (pCFJ90 92 ) at 5 ng/ μL, and hygromycin resistance (IR98 93 ) at 3 ng/μL. Independent lines were created from F1 animals selected for pharynx expression of the mCherry coinjection marker. Lines were maintained on Hygromycin selection plates as described above.
C-termin al GFP knock-in of his-72. C-terminal GFP knock-in of his-72 was performed as described previously using a self-excising selection cassette 89 .
Biosorter. Automated measurement of GFP negative animals in F1 and their F2 progeny. His-72::GFP was targeted with sg1, sg2, pool1 (sg2, 3, 4, 6, 8) or pool2 (sg3, 5,8). F1 generation was collected by bleaching 12 hrs after heatshock. These were either measured on the Biosorter fl ow system at larvae stage L3 or grown to adulthood to collect F2 generation which was then also measured at larvae stage L3.  35 , while the lethal vulva bursting occurs later after molting, in the adult stage 24 . RNA was chlorophorm-extracted as follows. Samples were thawed, 0.2 mL of chlorophorm added, incubated for 3 minutes, and centrifuged for 15 minutes at 12,000 x g at 4°C. The upper aqueous phase was transferred to a new tube, 2 μL GlycoBlue (30 μg) were added, 500 μL of isopropanol were added and sample was incubated for 10 minutes. Sample was centrifuged 10 minutes at 12,000 x g at 4°C, supernatant discarded, and 1 mL of 75% EtOH was added. Sample was centrifuged for 5 minutes at 7,500 x g at 4°C, supernatant removed, pellet air-dried and resuspended in 20 μL RNase-free water. RNA concentrations ranged between 1,000 -2,000 ng/μL, as determined on a Nanodrop ND-1000.

DNA sampling over generations (lin-41).
Mutated F1 samples were obtained as described above using large-scale mutagenesis by Cas9 heat shock induction. For this we used N2 as control and 3 lines with sgRNAs against the lin-41 3' UTR (sg15+sg16, sg26+sg27, sgPool). We conducted the experiment at 16°C and 24°C. 3000 L1 stage animals (F1 generation) were seeded on medium plates with OP50. After egg laying and hatching of the next generation (F2) after 3 or 5 days (24°C or 16°C) F1 and F2 were separated. For this, animals were washed from plates in a fi nal volume of 2 ml M9 buff er into 2 ml Eppendorf tubes. Adult animals sink faster and after circa 2-5 minutes are collected at the bottom of the tube, while L1 animals still swim. This was carefully monitored visually. When most adults (95%) had sunken to the bottom, supernatant M9, containing L1 stage animals, was removed to a separate tube. This was repeated three times by adding 2 ml M9 and separation by sinking. Adult animals were frozen for genomic DNA extraction in circa 20 uL M9. For generations F2-F4, 2000 L1 were seeded on new medium plates, and frozen as adults after separation from the next generation. Generation F5 was frozen at L1 stage. Genomic DNA extraction and targeted large amplicon sequencing was performed as described above.
Screen for regulatory sequences aff ecting phenotype. We targeted 8 genes with known RNAi-phenotypes (dpy-2, dpy-10, egl-30, rol-6, sqt-2, sqt-3, unc-26, unc-54) using diff erent sets of sgRNAs against regulatory regions. We used lines in which we targeted the 3′ UTR and for some genes we used additional lines targeting predicted enhancer, TATA-box, initiator (INR) and upstream/promoter regions. A list with all samples can be found in a supplementary table.
For each transgenic line we screened 35,000 -45,000 F2 animals produced from P0 with transiently induced Cas9 expression as described above. Animals were seeded onto NGM plates with food at a concentration of 15,000 per 15 cm plates or at 2,500 -5,000 per 10 cm plates. Plates were kept at 16°C or 24°C. We then directly screened these plates by eye. Additionally, we collected worms in M9 and dispensed worms in drops on an empty plate. We then observed worms moving in M9 and moving away after M9 was dried (<1 min.). Dpy, Unc, and Rol worms were identifi ed by morphology, their movement in M9 or slow and otherwise impaired movement away from the spot of dispension. Potential mutants were then picked and kept on plates for 2 -4 generations at 24°C to achieve homozygosity. Animals were then singled by phenotype and genotyped. This resulted in isolation of several mutant strains with the same genotype. We could not distinguish between cousins/ siblings coming from the same F1/F2 or independent mutants coming from independently mutated F1s. In these cases, we kept one representative strain. We determined that penetrance was complete for all alleles except for the sqt-2 locus (n>300 animals). For sqt-2 the penetrance varied between 10-100%. We scored the expressivity of the phenotypes into three categories (+, ++, +++) (n>300 animals). All the reported phenotypes have been determined and validated for several generations at 24°C. We also validated the absence of the extra-chromosomal transgenes judged by the red fl uorescent co-injection marker. For sqt-3 all isolated Dpy animals, characteristic for complete loss-of-function, contained large mutations aff ecting the coding frame. We therefore screened mainly for reduction-of-function alleles by screening for Rol animals. Non-Rol revertants of the sqt-3(ins) Rol animals were isolated using the small-scale approach on 6 cm plates with injection mixes imJJF215 or imJJF230.
Genotyping. Single worms were picked using a platin wire picking tool and immersed in 10 μL of worm lysis buff er (WLB) ( Pre-mRNA was specifi cally detected using intron-overlapping primers, while mRNA primers overlapped with exon-exon junctions. Controls without reverse transcriptase ("RT-") were done to ensure specifi c amplifi cation of cDNA and no amplifi cation from potential contaminating genomic DNA. Final values were obtained by normalizing to pre-mRNA or mRNA of tbb-2 and presented relative to N2 wild-type controls. Probes and primers can be found in a supplementary table. To prepare injection mixes, Cas9 protein was mixed with KCl and Hepes. crRNA and tracrRNA were annealed in duplex buff er by 5 min at 95 °C and ramp down to 25°C and added. Cas9/ tracRNA/crRNA mix was incubated at 37°C for 10 min. Then plasmids and ssDNA repair template were added and 10 P0 animals were injected. For each injection mix 8 F1s positive for the co-injection marker were picked and genotyped using two PCR reactions (one primer pair fl anking the insertion, the other with one primer binding in the insertion).

CRISPR-Downstream Analysis and Reporting Tool ("crispr-DART").
The source code along with installation instructions and sample input fi les can be accessed here: https://github.com/BIMSBbioinfo/crispr_DART. The source code for reproducing the fi gures generated based on the pipeline's output can be found here: https://github.com/BIMSBbioinfo/froehlich_uyar_et_al_2020.

Purpose of the crispr-DART pipeline.
In order to evaluate the outcomes of the CRISPR-Cas9 induced mutations by the protocol described in this study, we developed a computational pipeline to process the high-throughput sequencing reads coming from samples treated/untreated with CRISPR-Cas9. Although we developed the pipeline to address the hypotheses considering the specifi c experimental design in this study, we tried to make the pipeline as generic as possible to accommodate diff erent experimental setups, hoping that the pipeline can be useful to the scientifi c community carrying out genome editing experiments using the CRISPR-based technologies, in particular those that aim to introduce many combinations of mutations in a genome via inducing double-stranded DNA breaks repaired via the non-homologous end joining pathways. The pipeline can handle both short (single-or pairedend Illumina) and long reads (PacBio) from both DNA and RNA, so that the eff ects in the DNA editing can be also observed in the matching RNA samples. Each sample can contain multiple sgRNAs targeting multiple regions of the genome. A sample might contain reads coming from one or more individual genomes/transcriptomes. However, each sample is treated as if it is a single individual with a mosaic genome. The fi rst purpose of the pipeline is to serve as a quality control/reporting tool to evaluate the genomeediting experiment and address the following questions: Has the CRISPR-Cas9 treatment induced any mutations? If so, how are they distributed in the genome? Do the mutations that are commonly found in many reads originate at the intended cut site based on the designed guide RNA matching sites in the genome? How effi cient were diff erent guide designs in inducing DNA damage? Can we capture long deletions if there are multiple sgRNAs used in the same sample targeting nearby sites? How diverse are the deletions or insertions detected at the cut sites? We developed the pipeline to produce HTML reports collated into a website with interactive fi gures that help the user to quickly visualize and evaluate the outcomes of their experiment.
The second purpose of the pipeline is to produce many processed fi les containing information that can be useful for further analysis by external tools. Therefore, the pipeline's output consists of BAM fi les, bigwig fi les, BED fi les, and many diff erent tables containing information about insertions and deletions along with the reads in which they were detected. In this study, many of the fi gures made for the manuscript were generated based on these intermediate fi les to address the many custom questions.
Input description. crispr-DART is implemented using Snakemake 95 following the practices as implemented for the PiGx pipelines 96 . The input consists of a settings fi le in yaml format, which contains confi gurations for the tools used in the pipeline. Moreover, it contains fi le paths for where the sequencing reads are located, the target genome sequence to be used for mapping the reads, the sample sheet fi le (in comma-separated fi le format) which contains the experimental design, the fi le containing the genomic coordinates of the expected sgRNA cut sites (in BED fi le format), and a table (in tab-separated format) that is needed for when a pair of samples are to be compared (for instance to observe the diff erences in per-base distribution of deletions detected in a treated sample and an untreated control sample).
Steps of the pipeline. The pipeline consists of the following sequence of processing steps (see Supplementary Figure 2b): Pre-processing reads  Quality control (using fastqc 97 and multiqc 98 and quality improvement of reads (using Trim-Galore! 99 ) Mapping  Mapping/alignment of the reads to the genome (using BBMAP 100 ).
We use BBMAP for read alignment because it can handle both long and short reads, both single-end and paired-end reads, both DNA and RNA reads, and it can help detect both short and long insertion and deletion events.
 Re-alignment of reads with indels using GATK 101 . This step helps reconcile diff erent indel alignments to minimise the noise in alignments.

Finding indels and generating indel related statistics
 Extraction of indels from the BAM fi les using R packages    This table contains statistics about the effi ciency of each guide RNA in inducing mutations at the targeted site of the genome. The effi ciency of a sgRNA is defi ned as the ratio of the number of reads with an insertion/deletion that start or end at +/-5bp of the intended cut-site to the total number of reads aligned at this region.

Generating a website of HTML reports
All the pre-processed fi les from the previous steps are combined to generate interactive (where applicable) HTML reports from all the analysed samples that exist in the input sample sheet. For each targeted region (assuming a region of a few thousand base pairs that is sequenced), currently four diff erent reporting Rmarkdown scripts are run. The resulting HTML fi les are organised into a website using the `render_site` function of the Rmarkdown package 104 . Thus, all the processed data and outcomes can be quickly browsed through a website. The resulting website that is the output of the pipeline for this study can be browsed here: https://bimsbstatic.mdc-berlin.de/akalin/buyar/ froehlich_uyar_et_al_2020/reports/index.html How to reproduce the fi gures in the manuscript The source code and usage instructions for reproducing the fi gures in this manuscript, that were made based on the output of the crispr-DART pipeline, can be found here: https://github.com/BIMSBbioinfo/froehlich_uyar_et_ al_2020. Deletions were further categorized based on whether they overlap both let-7 micro-rna seed regions, and also those that don't overlap any of these defi ned regions.

RNA analysis: Comparing deletions in L1 and L4 stages
Deletion frequency values were computed and the ratio of deletion frequencies between L4 stage and L1 stage samples were computed in log2 scale. For each category of deletions, a wilcoxon rank-sum test was computed to test the null hypothesis that the stage specifi c abundance of deletions that overlap a let-7 binding site is not diff erent from those deletions that don't overlap any of these sites.

RNA analysis: Clustering pacbio reads and comparing genotype abundance in L1 and L4 stages
Reads from both L1 and L4 stage lin-41 RNA samples sequenced using PacBio that cover the region between chrI:9334840-9336100 (the region from the beginning of the amplifi ed segment up to the fi rst intron) were selected to make sure that all reads that go into analysis are covering the whole segment. For each read, the alignment of the read (including the inserted sequences) were obtained and all combinations of k-mers (k=5) were counted within these alignments allowing for up to 1 mismatch using Biostrings package 105 . Seurat package 106 was used to process the k-mer count matrix to do scaling, dimension reduction (PCA and UMAP) and networkbased spectral clustering. The clustering of long PacBio reads covering the region enabled us to cluster reads into genotypes, thus taking advantage of the length of the reads while also allowing for the high rate of indels in the PacBio reads (compared to the Illumina reads).

Fitness analysis
For this analysis, we have utilized lin-41 DNA samples sequenced (Illumina single-end sequencing) from multiple generations from F2 to F5 of the same pool of animals treated with sgRNA guides "sg15 and sg16", "sg26 and sg27" and "sg pool". We would like to address the question whether the deletions that exist at F2 are exposed to purifi ed selection over generations if they overlap the important sites on the 3'UTR region of lin-41. We do this analysis in two ways:  First, we count the deletions categorised by their overlap (or nonoverlap) with the important sites that exist in F2 generation and see how many of them still exist in later generations.
 Secondly, we do the same analysis at the level of reads: counting the reads with deletions that overlap or not overlap the important sites from generations F2 to F5. When comparing the number of reads, the read counts are normalized by the library sizes (total number of reads in the sample).
sgRNA effi ciency comparisons. For comparing observed effi ciencies to published prediction scores and other sgRNA characteristics, these were manually extracted from the CRISPOR web application (http://crispor.tefor. net/) 107 . Browser shots. Browser shots were compiled using indel profi les and top indels provided by the computational pipeline as BigWig and BED fi les and loading them into the UCSC genome browser 108 or the IGV browser 109 followed by export as vector graphics compatible format. We used C. elegans genome version ce11/WBcel235 including 26 species base-wise conservation (PhyloP).
Graphics. Final fi gures were assembled using Adobe Illustrator without changing the plotted data, with adjustments to composition, axis labels, line widths, colors, etc.
Manuscript formatting. This manuscript was composed using Adobe InDesign starting with a template from https://github.com/cleterrier/ ManuscriptTools. We thank Christophe Leterrier for sharing it.

Data reporting.
No statistical methods were used to predetermine sample size. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.
Strain and plasmid availability. C. elegans mutant strains will be available through the Caenorhabditis Genetics Center (CGC) as indicated in a supplementary table. Plasmids generated for this work for heat shock Cas9 expression (pJJF152) and proof of concept sgRNAs (targeting SECGFP, dpy-10, sqt-3) will be available from Addgene as indicated in a supplementary table or for specifi c sgRNAs upon direct request from the authors.