Pervasive epistasis in cell proliferation pathways modulates neurodevelopmental defects of autism-associated 16p11.2 deletion

Rare CNVs such as the 16p11.2 deletion are associated with extensive phenotypic heterogeneity, complicating disease gene discovery and functional evaluation. We used RNA interference in Drosophila melanogaster to evaluate the phenotype, function, and interactions of conserved 16p11.2 homologs in a tissue-specific manner. Using a series of quantitative methods for assessing developmental and neuronal phenotypes, we identified multiple homologs that were sensitive to dosage and showed defects in cell proliferation, including KCTD13, MAPK3, and PPP4C. Leveraging the Drosophila eye for studying gene interactions, we performed 565 pairwise knockdowns of gene expression, and identified 24 interactions between 16p11.2 homologs (such as MAPK3 and CORO1A, and KCTD13 and ALDOA) and 46 interactions with other neurodevelopmental genes (such as MAPK3 and PTEN, and KCTD13 and RAF1) that significantly enhanced or suppressed cell proliferation phenotypes. Integration of fly interaction and transcriptome data into a human brain-specific genetic network allowed us to identify 982 novel interactions of 16p11.2 homologs, which were significantly enriched for cell proliferation genes (p=3.14×10−12). Overall, these results point towards a new model for pathogenicity of rare CNVs, where CNV genes interact with each other in conserved pathways to modulate expression of the neurodevelopmental phenotype.


INTRODUCTION
Rare recurrent copy-number variants with breakpoints typically mapping within segmental duplications are a significant cause of neurodevelopmental disorders, such as intellectual disability/developmental delay (ID/DD), autism, epilepsy, and schizophrenia 1 . Gene discovery within rare syndromic CNVs has traditionally involved mapping the disease-associated region using atypical CNVs, inversions, or translocations to identify a causative gene that explains the distinct phenotypes associated with the CNV, followed by detailed functional evaluation of that gene using animal models. Using this approach, the retinoic acid induced 1 gene (RAI1) was identified as the locus responsible for the core features of Smith-Magenis syndrome 2 , and individual genes within chromosome 7q11.23 were connected to specific Williams-Beuren syndrome phenotypes, such as ELN for cardiovascular features 3 . Absence of atypical deletions for other CNVs required more direct functional evidence for implicating a candidate gene. For example, the role of TBX1 in aortic arch defects observed in individuals with 22q11.2 deletion/DiGeorge syndrome was identified through a functional screen for cardiac features in a series of mouse models carrying overlapping deletions of the human syntenic region 4 . All of these examples provide evidence that dosage alteration of one or more genes within the syndromic CNV interval contribute to the observed phenotypes.
Unlike rare CNVs associated with a consistent set of phenotypes, more recently described rare CNVs are associated with a range of neurodevelopmental features, and are also reported in unaffected or mildly affected individuals 1 . One such CNV is the 16p11.2 deletion, which encompasses 593 kbp and 25 unique genes. The deletion was originally identified in individuals with autism 5,6 , and subsequently reported in children with ID/DD 7 , epilepsy 8 , and obesity 9 .
Several themes have emerged from recent studies on dissecting the role of individual genes within the 16p11.2 deletion towards neurodevelopmental phenotypes. First, extensive heterogeneity and incomplete penetrance of the associated phenotypes adds additional challenges to genetic mapping strategies that use atypical variants. Second, while this deletion is enriched within various neurodevelopmental disease cohorts, exome sequencing studies of hundreds of individuals have not identified any individual genes within this region as causative for these diseases on their own [10][11][12][13] . Third, functional studies using cellular 14,15 , mouse [16][17][18][19] , and zebrafish models [20][21][22] have implicated several different genes within 16p11.2 in neurodevelopmental phenotypes. These findings suggest that the observed phenotypes in 16p11.2 deletion are not caused by haploinsufficiency of a single causative gene, but rather are modulated by multiple dosage-sensitive genes in the region, potentially through combinatorial mechanisms within pathways related to neurodevelopment. This model is also consistent with a recent observation that pathogenic CNVs are more likely to contain clusters of functionally related genes than benign CNVs 23 , suggesting intra-CNV genetic interactions as a potential cause for CNV pathogenicity. Therefore, an approach that combines a systematic functional evaluation of each gene within 16p11.2 and its genetic interactions is necessary to identify key neurodevelopmental pathways and molecular mechanisms of disease.
Evaluation of gene interactions in neurodevelopment requires a system that is sensitive to genetic perturbations but, at the same time, allows for performing interaction studies in the nervous system without compromising viability of the organism. Drosophila melanogaster provides such a model, as developmental processes, synaptic mechanisms, and neural structure and signaling are conserved between flies and vertebrates 24 . In fact, several neurodevelopmental disorders, such as Angelman Syndrome 25 , Rett Syndrome 26 , and Fragile X syndrome 27 , have been modeled in flies. We used the power of Drosophila melanogaster as a genetic model to perform a series of quantitative and high-throughput assays to systematically characterize phenotypes, function, cellular mechanisms, and interactions of conserved homologs of human 16p11.2 genes. Our data suggest a pervasive epistatic model for disease pathogenicity, where multiple 16p11.2 genes are sensitive to dosage imbalance and participate in complex interactions that enhance or suppress the phenotypic effects of each other within cellular proliferation pathways, and in turn are modulated by other genes in the genetic background.

Multiple 16p11.2 homologs contribute to neurodevelopmental phenotypes
We identified 14 fly homologs from the 25 human 16p11.2 genes (Table S1), and used 31 RNA interference (RNAi) lines and tissue-specific GAL4 drivers to knockdown the expression levels of individual homologs ubiquitously or in neuronal, eye, or wing tissues (Figures 1 and 2A).
RNAi is an effective strategy to model partial reduction of gene expression, which in principle recapitulates the effect of a heterozygous microdeletion, and for high-throughput screening of genes for tissue specific phenotypes. We used multiple independent UAS-RNAi transgenes targeting the same gene to validate our results ( Figure S1A, Table S2), and used stringent quality control to eliminate lines that showed phenotypes due to off-target or positional effects ( Figure   S1A). Using quantitative PCR, we measured the reduction in gene expression for each line with neuronal knockdown (using Elav-GAL4>Dicer2 at 25°C), and on average achieved approximately 50% reduction in gene expression for the tested 16p11.2 homologs ( Figure S1B).
We performed a series of quantitative assays on 16p11.2 homologs for more than 20 phenotypes that have been classically used to measure conserved developmental function in flies, and identified lethality and a variety of morphological phenotypes due to ubiquitous and tissuespecific knockdown 27 (Figures 1 and 2A). For example, seven homologs were lethal at the larval or pupal stage with ubiquitous knockdown, indicating that these genes are essential for viability and development in Drosophila 28 . We next performed pan-neuronal knockdown experiments and tested for several nervous system phenotypes, such as climbing assays for motor impairment, spontaneous seizures, and head size. We performed negative geotaxis experiments to measure locomotor function and identified dramatic reductions in the climbing ability of ALDOA ald and MAPK3 rl knockdown flies throughout the testing period ( Figure 2B, Figure S2A). Since about 20% of individuals with 16p11.2 deletion manifest seizures 29 , we next used a recently developed spontaneous seizure assay that assesses unprovoked seizures in their native state, which better recapitulates human seizures in Drosophila 30 . We found that MAPK3 rl , CDIPT pis , PPP4C pp4-19C , and KCTD13 CG10465 knockdown flies were more likely to show seizure phenotypes compared to controls ( Figure 2C, Figure S2B). We further measured fly head sizes, as cellular signaling pathways such as Hippo and TOR signaling contribute to brain development 31 and individuals with the deletion show macrocephaly features 29 . Using CTNNB1 arm and PTEN dpten knockdown flies as positive controls 32 , significant changes in head sizes were found in certain 16p11.2 homologs, including an increased head size for KCTD13 CG10465 knockdown similar to that reported in zebrafish models 21 as well as increased head sizes for C16ORF53 pa1 and CDIPT pis ( Figure S2C).
We further examined deeper cellular features, including neuromuscular junction, dendritic arborization, and axonal targeting, to understand the molecular basis of the observed neuronal features. Drosophila neuromuscular junction (NMJ) is a well-established model for studying synapse growth defects, and alterations in NMJ architecture have been documented in genes associated with autism 27 . We found significant differences in NMJ structure with knockdown of CDIPT pis and FAM57B CG17481 compared to controls, suggesting altered growth and development of the NMJ in these flies ( Figures 2D and S2D). The architecture of dendritic arbors also plays an important role in neural circuit formation, and defects in dendrites are associated with neurodevelopmental disorders such as schizophrenia and autism 33 . To assay dendritic growth and structure, we examined large branched dendrites of the class IV ddaC sensory neuron in intact larvae after gene knockdown with the ppk-GAL4 driver 34 , and observed decreased complexity in dendritic arborization with knockdown of KCTD13 CG10465 and TAOK2 dTao-1 ( Figures 2E, Figures S2E and S2F). Another hallmark of nervous system development is the accuracy of synaptic connections, which is determined by the guidance of axons to their correct targets 35 . We explored axonal targeting by staining larval eye discs of flies using chaoptin antibody, and observed aberrant targeting in KCTD13 CG10465 and MAPK3 rl flies ( Figure 2F). In summary, we found multiple developmental and neuronal defects for each of the homologs, indicating the pleiotropic effect of conserved 16p11.2 genes and their importance in neurodevelopment.

Drosophila eye models identify cellular proliferation phenotypes in 16p11.2 homologs
Decades of studies have shown that the Drosophila eye is an accessible and sensitized experimental system for quantitative studies of nervous system development and function, as genetic defects that alter the development of a single cell type can lead to observable rough eye phenotypes 36,37 (Figures 3A and 3B). In fact, genetic interaction studies using the fly eye have led to the discovery of novel modifier genes in nervous system disorders, such as Rett syndrome and spinocerebellar ataxia type 3 26,38 , as well as conserved developmental processes 39 . To quantify the degree of severity of the eye phenotypes, we developed and tested a computational method called Flynotyper that calculates a phenotypic score based on the disorderliness of ommatidial arrangement at high sensitivity and specificity 40 . We performed eye-specific knockdown of gene expression using the GMR-GAL4 driver with and without Dicer2 for all Drosophila 16p11.2 homologs, and compared the degree of phenotypic severity as measured by Flynotyper to controls with the same genetic background ( Figures 3C and 3D, Figure S3A-C).
We found a strong correlation (Pearson correlation, r=0.69, p=2.88×10 -6 ) between the percentile ranks of all tested RNAi lines with Dicer2 and without Dicer2 ( Figure S3D). As shown in Figure   3D, we observed a range of severe but significant eye defects for nine homologs, which were comparable to that of genes associated with neurodevelopmental disorders such as CHD8 kis , SHANK3 prosap , SCN1A para , and PTEN dpten (Table S3, Figure S3D). For example, the severity of eye defects in KCTD13 CG10465 , DOC2A rph , and PPP4C pp4-19C knockdown flies had phenotypic scores greater than the 85 th percentile of the tested 39 fly homologs of human neurodevelopmental genes (Table S3). These results suggested that knockdown of 16p11.2 homologs affect development of the fly eye to varying degrees of severity, which mirrors the global developmental and neurological defects observed with ubiquitous and pan-neuronal knockdown.
To investigate the cellular basis of the rough eye phenotype observed in individual gene knockdowns, we stained the pupal eye imaginal disc with Discs large (Dlg) antibody for ommatidial cells and Phalloidin for photoreceptor neurons, and screened for anomalies in different cells in the developing eye ( Figures 3A-3C). A variety of cellular defects leading to altered structure of the hexagonal lattice were observed with knockdown of seven of the homologs, suggesting potential alterations in cellular proliferation (Figures 3E and 3F, Figure   S4). For example, KCTD13 CG10465 knockdown flies showed a drastic increase in the number of cone cells, secondary pigment cells, and photoreceptor neurons, while MAPK3 rl knockdown showed a decreased number of interommatidial cells and photoreceptor neurons, with a consequent loss of the hexagonal structure in the ommatidia. Similarly, ALDOA ald knockdown flies had misplaced bristle cells as well as an increase in secondary and photoreceptor cells, while PPP4C pp4-19C knockdowns showed severe rotational defects and a complete loss of the ommatidial architecture. Overall, we found that knockdown of several 16p11.2 homologs contribute to defects in cell count and patterning of different cell types, including photoreceptor neurons and ommatidial cells, during development.
We further investigated cellular mechanisms associated with the observed developmental defects, as recent functional studies have implicated defects in neuronal proliferation as a cellular mechanism underlying autism disorders 41 . In fact, genome-wide CNV and exome sequencing studies of individuals with autism have uncovered pathogenic variants enriched for cell proliferation genes 42,43 . We used phospho-Histone H3 (pH3) antibody and bromodeoxyuridine (BrdU) staining to identify dividing cells, and counted the number of stained cells posterior to the morphogenetic furrow of the developing larval eye ( Figures 3A and 3C, Figure S4). Several homologs showed a significant alteration in dividing cell counts ( Figure 3G). For example, we found an increase in mitotic cell count with knockdown of KCTD13 CG1046 , CDIPT pis , and ALDOA ald , while MAPK3 rl knockdown flies showed a significant reduction in proliferating cells.
No changes in cell differentiation were observed using Elav staining in KCTD13 CG10465 and MAPK3 rl knockdowns, suggesting that these two genes are specifically involved in cell proliferation (Figures S4B). Consistent with the proliferation phenotypes, we also observed an overall increase in the adult eye area in four RNAi knockdowns comparable to flies with knockdown of PTEN dpten , a known cell proliferation gene 44 , as well as a decrease in eye area for MAPK3 rl flies ( Figure 3H). Overall, our analysis of individual gene knockdowns showed that reduced expression of individual 16p11.2 homologs cause defects in cell proliferation and organization.

Interactions between 16p11.2 homologs modulate neurodevelopmental and cellular phenotypes
Our phenotypic and functional studies of individual 16p11.2 homologs showed that many genes within the CNV region are involved in neurodevelopment, indicating that no single gene in the region is solely responsible for the observed neuronal phenotypes of the deletion. Based on these observations we hypothesized that interactions between genes within the 16p11.2 region may contribute to the observed phenotypes. To systematically test interactions between 16p11.2 homologs in the developing fly eye, we selected a subset of four homologs, including PPP4C pp4-19C , MAPK3 rl , KCTD13 CG10465 and DOC2A rph , and generated recombinant lines expressing their respective RNAi lines with the GMR-GAL4 driver. We selected these genes as primary drivers of neurodevelopmental phenotypes based on severity of phenotypes with various tissue-specific knockdowns (Figure 2A), published functional studies in mouse and zebrafish (Table S1), and identifiable eye phenotypes amenable for large-scale modifier screens ( Figure 3D).
We generated 52 two-locus fly models (123 total lines) by reducing gene expression of each of the four homologs in combination with the 13 other 16p11.2 homologs. We used manual eye scoring and Flynotyper to compare these pairwise knockdown lines to respective control flies with single gene knockdowns. In this way, we identified 24 (Figures 4B and S5B). Further, the glossy eye phenotype observed in one-hit knockdown of PPP4C pp4-19C was suppressed by reduced expression of five homologs, including YPEL3 CG15309 and CDIPT pis ( Figures 4C and S5C). We also observed an enhancement of the rough eye phenotype with double knockdown of PPP4C pp4-19C and KCTD13 CG10465 at 30°C, which we initially suspected to be due to the severity of KCTD13 CG10465 knockdown by itself. To dissect this interaction, we performed reciprocal crosses of KCTD13 CG10465 RNAi lines with PPP4C pp4-19C at 25°C, and confirmed the phenotypic enhancement observed at 30°C ( Figure   S5D). Finally, the rough eye phenotypes in DOC2A rph knockdown models were suppressed by reduced expression of CDIPT pis and ALDOA ald ( Figure S5E).
We tested a subset of 16p11.2 interaction pairs for alteration in their cellular phenotypes by staining the pupal and larval eye discs with anti-Dlg and anti-pH3, respectively ( Figures 4E   and S5F). Assessing the count, structure, and orientation of the cells in the developing eye discs confirmed several interactions documented in the adult eyes. For example, reduced expression of ALDOA ald in MAPK3 rl models rescued the rotation errors and primary cell defects in the pupal eye, as well as proliferation defects in the larval eye ( Figures 4E-4G, Figure S5H). Similarly, ALDOA ald knockdown suppressed the cone cell defects, secondary cell defects, and rotation errors observed in KCTD13 CG10465 knockdown pupae, and showed a significant reduction in the number of proliferating cells compared to the KCTD13 CG10465 larval eye ( Figures 4E-4G, Figure   S5H). Although reduced expression of TAOK2 dTao-1 did not rescue external eye defects in MAPK3 rl and KCTD13 CG10465 knockdown flies, we observed partial rescue of cellular defects in the pupa and a significant rescue of proliferation defects in MAPK3 rl and KCTD13 CG10465 larval eyes ( Figures 4E-4G, Figure S5F-S5H). To test if the two-hit interactions observed in the fly eye were also relevant in the nervous system, we evaluated the accuracy of retinal axon innervation into the lamina and medulla of the brain using anti-chaoptin in larvae with knockdown of both KCTD13 CG10465 and CORO1A coro , and confirmed a complete rescue of the axonal targeting defects observed in single-hit KCTD13 CG1046 flies ( Figure 4D). Some literature evidence exists for the functional interactions documented in our study.
For example, MAPK3 and TAOK2 are both involved in synaptic assembly and signaling 45 , and ALDOA was identified as a member of the MAP kinase (ERK1/2) interactome in differentiating epidermal and neuronal cells 46 . We also found that the tested 16p11.  Table S5). Interestingly, 17 of these interactions of 16p11.2 homologs were with genes known to be involved in cell proliferation. For example, knockdown of the Wnt signaling pathway gene 48 CHD8 kis resulted in a complete rescue of MAPK3 rl phenotype as well as a strong suppression of the KCTD13 CG10465 rough eye phenotype ( Figure   5C). Similarly, reduced expression of beta catenin, CTNNB1 arm , significantly enhanced the phenotypes observed with MAPK3 rl , KCTD13 CG10465 , and PPP4C pp4-19C one-hit knockdowns ( Figure 5C). Knockdown of SHANK3 prosap , a gene that codes for a post-synaptic scaffolding protein and is associated with autism 49 , suppressed the rough-eye phenotype of KCTD13 CG10465 , MAPK3 rl , and PPP4C pp4-19C flies ( Figure 5C). These data show that multiple 16p11.2 homologs interact through conserved neurodevelopmental genes that potentially act as intermediates, whose knockdown modulates the expression of the ultimate phenotype. Interestingly, six genes within the 16p11.2 distal CNV region 50 interacted with 16p11.2 homologs. For example, reduced expression of SH2B1 lnk fully rescued the rough-eye phenotype observed with knockdown of MAPK3 rl , while ATXN2L atx2 knockdown led to a more severe phenotype in PPP4C pp4-19C , KCTD13 CG10465 , and DOC2A rph knockdown flies ( Figure 5C). These results suggest overlapping functional roles in neurodevelopment for genes within the proximal and distal 16p11.2 regions apart from chromosomal contacts between the two syntenic segments in humans 51 . We further assessed the cellular mechanisms responsible for suppression of the MAPK3 rl rough eye phenotype by simultaneous knockdown of the autism-associated tumor suppressor 52 PTEN dpten ( Figures 5C-5F). We observed a complete rescue of the bristle group defects, rotation errors, primary and secondary cell defects, and photoreceptor cell counts in the pupal eye of the two-hit knockdowns, and a complete rescue of cell proliferation defects observed with MAPK3 rl singlehit knockdown (Figures 5D-5F, Figure S6E and 6F).

Transcriptome studies of 16p11.2 homologs show a network of genes enriched for cell proliferation
To examine functional processes associated with 16p11.2 homologs, we selected six 16p11.2 homologs, MAPK3 rl , KCTD13 CG10465 , DOC2A rph , CORO1A coro , C16ORF53 pa1 , and CDIPT pis , based on high phenotypic severity, and performed RNA sequencing on fly heads for each knockdown to identify differentially-expressed genes (Table S6). We first conducted parametric gene-set enrichment analysis 53 to identify Gene Ontology terms enriched for human homologs of up-or down-regulated genes in each fly model relative to the control ( Figure S7A, Table S6).
Several terms related to neurodevelopment, including cell cycle process, neurogenesis, neuron differentiation, and cell-cell signaling, were significantly enriched in each knockdown model (p<0.01, corrected by Benjamini-Hochberg method) ( Figure S7A). We also observed a significant enrichment for intellectual disability genes (p=0.042 after Bonferroni correction, onetailed Student's t-test) among genes differentially expressed by KCTD13 CG10465 knockdown ( Figure S7B). Based on the cell proliferation phenotypes observed in the fly eye, we further constructed a network of differentially-expressed genes annotated for cell cycle and proliferation in humans ( Figure S7C). Interestingly, we found a significantly high degree of overlap among differentially-expressed cell proliferation genes between the knockdown models (empirical p<0.001), with 37.2% (64/172) of these genes differentially expressed in two or more models and 20.3% (35/172) differentially expressed in at least three knockdown models. These results provide additional evidence that the 16p11.2 homologs function in well-connected cell proliferation pathways in both Drosophila and humans.
We next selected 22 genes that were among the most up-regulated genes in KCTD13 CG10465 and MAPK3 rl models for two-hit interaction experiments, and identified 18 genes whose knockdown suppressed the rough eye phenotypes in the MAPK3 rl and KCTD13 CG10465 flies ( Figures 6A and 6B). For example, knockdown of COX6AL coxA2 fully rescued the MAPK3 rl rough eye phenotype, and knockdown of CNGA2 CG42260 , CYP24A1 cyp12d1-d and RAF1 CG14607 suppressed the KCTD13 CG10465 phenotype. We further examined the cellular phenotypes of  Figure 6F). Of note, although the transcriptome analysis was performed on fly brains using a neuron-specific driver, we were able to validate those interactions in the fly eye, supporting the utility and veracity of using the Drosophila eye to study nervous system interactions.

Concordance of functional interactions with human brain-specific genetic interaction network provides translational relevance
To assess the relevance of the identified functional interactions in our fly experiments to neurodevelopment in humans, we explored the functional context of the human 16p11.2 genes and their involvement in cell proliferation, specifically in relation to brain biology. Using a Bayesian network of predicted and known genetic interactions in the brain, we first mapped 14 homologs of 16p11.2 genes and 35 interacting genes identified from fly experiments onto a human brain-specific gene interaction network 54,55 , and then identified additional genes in the network that connected these genes to each other. Overall, we found 982 brain-specific interactions, with 39 out of the 49 tested genes connected through 428 novel genes within the network ( Figure 6G, Table S7). A significant enrichment for cell cycle and cell proliferation function was identified among these novel connector genes (96/428, one-tailed Fisher's exact test, p=3.14×10 -12 ). Additionally, the connector genes were enriched for genes related to neurodevelopment, including FMRP-protein binding genes (p=3.34×10 -14 , one-tailed Student's ttest) and genes involved in post-synaptic density (p=3.31×10 -32 ), as well as genes differentially expressed in the fly knockdown models (p=0.0215). These results suggest a strong concordance between data obtained from fly two-locus experiments with putative interactions identified in the human nervous system, and provide a novel set of candidates that could be potential therapeutic targets for the deletion phenotypes.

DISCUSSION
We used the sensitive genetic system of Drosophila to identify conserved functions and interactions of 16p11.2 homologs that are relevant to the human 16p11.2 deletion. First, although we did not test every 16p11.2 gene, we identified multiple interactions of conserved 16p11.2 homologs that are consistent with published biochemical studies as well as functional gene networks constructed from human co-expression and protein-protein interaction datasets.
Second, the composite system of the fly eye allowed us to assay multiple genes and hundreds of interactions using high-throughput and quantitative assays for neurodevelopment without compromising the organism's viability. In fact, we were able to validate nervous system-specific interactions identified through transcriptome studies using both the fly eye and a human brainspecific network, which provides strong support for the use of genetic screens in the Drosophila eye for studying mechanisms of disease in the nervous system. Our results provide evidence for specific interactions in this region that can be investigated further in more sophisticated neurobiological systems, such as human stem cells, mouse, and zebrafish. Third, we observed conserved cellular mechanisms for pathogenicity of neurodevelopmental phenotypes in Drosophila, such as cell proliferation. These results point towards cellular proliferation pathways as a potential therapeutic target for 16p11.2 deletion phenotypes and related neurodevelopmental disorders. For example, rapamycin has been shown to rescue cellular and behavioral phenotypes in mouse knockouts of PTEN 56 . In the context of the deletion, the identified suppressors of multiple 16p11.2 homologs both within and outside the region, such as ALDOA, CORO1A, CHD8, PTEN, and RAF1, could be targeted for therapeutic strategies. Overall, our data show that Drosophila models of neurodevelopment can provide a new dimension to identifying gene interactions and cellular mechanisms of pathogenicity, and can be integrated with data from other animal and cellular models to fully explain the complex basis of neurodevelopmental disorders in humans.
Our results suggest that multiple 16p11.2 homologs contribute to a range of phenotypes that have key roles in different tissue types and organ systems, indicating pleiotropic effects in Drosophila that mirror the multitude of phenotypes observed in humans. These data are consistent with other functional studies of 16p11.2 knockdown or knockout models in mouse and zebrafish, which show abnormal neuronal or developmental phenotypes (Table S1).
Additionally, several 16p11.2 genes have similar intolerance to variation and likelihoods of having loss-of-function mutations compared to causative genes in syndromic CNVs 1 , including RAI1, SHANK3, and NSD1. The presence of multiple genes in 16p11.2 that are individually potentially pathogenic and behave similarly to classical neurodevelopmental genes, as determined by RVIS 57 and pLI scores 58 , suggest that interactions between these genes are necessary to modulate the effects of each gene in the deletion. Further, the additive effects of haploinsufficiency of all 16p11.2 genes cannot completely explain the clinical features of the deletion, as all patients with the deletion should then manifest some degree of the affected phenotype. However, this is not the case, as several clinical features of 16p11.2 deletion are not completely penetrant, including autism 59 and macrocephaly 29 . Based on our results, we propose a novel epistatic model for the pathogenicity of 16p11.2 deletion, where the phenotypic effects of the whole deletion may not equal the sum of the phenotypic effects due to disruption of its constituent genes. Rather, the interactions between the genes within the deletion, acting through common cell proliferation pathways, determine the phenotypic severity ( Figure 7). These interactions can suppress or rescue the one-hit phenotypes, providing evidence for epistasis as a phenomenon for the incomplete penetrance of clinical features associated with the deletion. The phenotypic variability of the deletion can therefore be explained by variants in upstream regulatory regions or modifier genes elsewhere in the genome that also participate in common pathways ( Figure 7). This epistatic model is in contrast to that reported for syndromic CNVs, where the core phenotypes can be due to a single gene (such as RAI1 in Smith-Magenis syndrome) or a subset of individual genes in the contiguous region (as in Williams syndrome).
While this model is likely applicable to other variably-expressive CNVs, our results further suggest the importance of genetic interactions towards causation and modulation of neurodevelopmental disease. Our study emphasizes the need for a function-based analysis in addition to sequencing studies towards discovery of gene function in the context of genetic interactions.

Identification of Drosophila homologs
We first queried the Drosophila genome for homologs using DRSC Integrative Ortholog Prediction Tool (DIOPT) 60 , reciprocal BLAST and ENSEMBL database 61 for each of the 25 human 16p11.2 genes (Table S1). We narrowed down the list of homologs to 15 genes with a DIOPT score of 3 or greater, or in the case of KIF22 klp68D , high query coverage and percentage identity in BLAST. Of the 15 selected homologs, RNAi lines were available for all homologs except INO80E, and therefore we were unable to characterize this gene. Therefore, 14 homologs of 16p11.2 genes were used in this study. We used a similar strategy for identifying homologs for other genes tested for interactions in this study. As this study is focused on studying the functional role of human genes in a fly model, we represent the identified fly homologs in the format of HumanGene FlyGene (for example, MAPK3 rl ). We confirmed that all tested 16p11.2 homologs were expressed in the fly eye using database and literature searches [62][63][64][65] .

Drosophila stocks and genetics
Tissue-specific knockdown of 16p11.2 homologs and other genes tested in this study was achieved with the UAS-GAL4 system 66 Table S2.
To study two-locus models, we first generated individual fly stocks with reduced expression for KCTD13 CG10465 , MAPK3 rl , PPP4C pp4-19C , and DOC2A rph containing GMR-GAL4 and UAS-RNAi on the same chromosome. We tested to ensure there was adequate GAL4 to bind to two independent UAS constructs in the two-locus models ( Figure S1E). Females from the stocks with constitutively reduced gene expression for each of these four genes were then crossed with other RNAi lines to achieve simultaneous knockdown of two genes in the eye.  (Tables S4 and S5).

RNA extraction and quantitative real-time PCR
We assessed mRNA expression by performing quantitative real-time PCR (qRT-PCR) experiments on cDNA samples isolated from fly heads. Briefly, RNAi lines were crossed with Elav-GAL4 driver at 25°C, and F1 female progeny were collected in groups of 40-50, quickly frozen in liquid nitrogen, and stored at −80°C. For RNA extraction, the heads were separated from their bodies by repetitive cycles of freezing in liquid nitrogen and vortexing. Total RNA was isolated using TRIZOL (Invitrogen), and reverse transcription was performed using qScript cDNA synthesis kit (Quanta Biosciences). RNA was also isolated from fly heads from GMR-GAL4 crosses for a subset of genes to compare the gene expression with fly heads from Elav-GAL4 crosses ( Figure S1C). Quantitative RT-PCR was performed using an Applied Biosystems Fast 7500 system with SYBR Green PCR master mix (Quanta Biosciences) using optimized protocols. A list of primers used for the qRT-PCR experiments is provided in Table S8.

Quality control
We checked the insertion site of the RNAi constructs to identify and remove any fly lines that may show phenotypes due to insertion-site effects ( Figure S1A). While RNAi transgenes for the Bloomington lines are inserted at the attP2 site on chromosome 3L with no expression or effect on the nervous system, thorough analysis of RNAi lines obtained from VDRC stock center was required to rule out off-target effects. We obtained 2 types of lines from VDRC: GD lines, which are P-element-based transgenes with random insertion sites, and KK lines, which are phiC31-based transgenes with defined insertion sites 28 . In order to rule out any effect of insertion of the RNAi construct in the GD lines, we mapped the insertion site by performing Thermal Asymmetric Interlaced PCR (TAIL-PCR) and Sanger sequencing. The TAIL-PCR method was modified from a protocol developed in B. Dickson's lab, based on published protocol 68 . The first round of PCR was performed with a 1:100 dilution of a genomic DNA preparation with Taq polymerase using three degenerate forward primers (AD1, AD2 and AD3) and a specific reverse primer (T1BUAS) (see Table S8 for TAIL-PCR primers). The second PCR reaction was set up using 1:50 dilution of the first PCR as template, with the AD primer as the forward primer and T2D as the specific reverse primer. The second PCR products were then visualized in 1%  69 . We observed non-specific shriveled wings in three out of seven KK lines of 16p11.2 homologs with Elav-GAL4, and these three lines also showed increased expression of tiptop ( Figure S1D). Therefore, we excluded these KK lines from neuronal experiments using Elav-GAL4. However, we found that overexpression of tiptop (using UAS-tio) with GMR-GAL4 showed a rough eye phenotype and reduced pigmentation confined to the right side of the eye, distinct from the eye phenotypes observed in the KK lines ( Figure S1D).
Further, we did not observe any changes in the expression of numb in the fly lines used in this study ( Figure S1D).

Climbing assay
Fly crosses were set up at 25°C with Elav-GAL4 to achieve neuronal knockdown. Four genes, PPP4C pp4-19C , ALDOA ald , TAOK2 dTao , and KIF22 klp68D , showed lethality when neuronal expression was reduced using RNAi at 25°C, and therefore were tested at room temperature.
KIF22 klp68D lines were also lethal when raised at room temperature. For each genotype, groups of 10 flies were transferred to a climbing vial and tapped down to the bottom. They were allowed to climb past a line marked at 8 cm from the bottom of the vial, and the number of flies crossing the 8 cm mark at 10 seconds was recorded. For each group, this assay was repeated nine more times with one-minute rest between each trial. These sets of 10 trials for each group were repeated daily for ten days, capturing data from flies aged day 1 to day 10. All experiments were performed during the same time of the day for consistency of the results. Two-way ANOVA and pairwise two-tailed T tests were used to determine significance for each genotype and day of experiment (Table S9).

Spontaneous seizures assay
Newly eclosed flies of the relevant genotypes were collected and aged for 7 days. Male and female flies were isolated at least 1 day after collection to ensure all females had mated. After aging, flies were transferred individually into the chambers of a 4×5 mating plate using a manual aspirator. The plate was then placed on a standard light box, where the flies were allowed to acclimate for 5 min. Fly behavior was recorded at 30 frames/sec for 5 min using a Canon High Definition Vixia HFM31 Camcorder (resolution 1920 x 1080). Each fly's behavior during the viewing window was then assessed for abrupt, involuntary seizure-associated movements, which manifest as rapid repositioning of the flies within the chamber as previously described 30 . For each knockdown genotype, the total number of flies that exhibited spontaneous seizure events was assessed and compared to control flies using a Fisher's exact test. We next determined the number of seizing events per seizing fly, comparing knockdowns with controls using a Student's t-test (Table S9).

Dendritic arborization assays
RNAi lines were crossed to a UAS-Dicer2; ppk-GAL4, UAS-mCD8-GFP driver at 25ºC, and embryos were collected at 24 hours on apple juice plates. First instar larvae, eclosed from the embryo, were transferred to the food plate and allowed to age for 48 hours at 25ºC before live imaging. Third instar larvae were collected, washed in PBS, and transferred dorsal side up to a glass slide containing a dried agarose pad with a coverslip on top secured with sticky tape. Zstack images of Class IV Dendritic Arborization neurons were acquired using a Zeiss LSM 800 confocal microscope and processed using ImageJ 70 to a scale of 5.0487 pixels/micron. Using an in-house Java plug-in, four concentric circles with a distance of 25 microns between each circle were placed on the images, with the cell body as the center. A manual Sholl analysis was conducted by counting the number of intersections of dendritic branches on each of the concentric circles. Total and average number of intersections were calculated and normalized to the width of the hemisegment of each sampled neuron to control for slight variation in larval sizes. Two-way ANOVA and pairwise two-tailed T tests were used to determine significance of the number of intersections in each genotype and concentric circle, and two-tailed Mann-Whitney tests were used to determine significance of the total number of intersections (Table   S9).

Phenotypic analysis of fly eyes using Flynotyper
We used GMR-GAL4 drivers with and without Dicer2 to achieve eye-specific knockdown, and imaged 2-3 day old flies using an Olympus BX53 compound microscope with an LMPlanFL N 10X 0.25 NA air objective (Olympus, Tokyo, Japan), at 0.5X magnification and a z-step size of 12.1µm. We used CellSens Dimension software (Olympus Optical) to capture the images, and stacked the image slices using Zerene Stacker (Zerene Systems, USA). All eye images presented in the figures are maximum projections of consecutive 20 optical z-sections. Eye area was calculated from each image using ImageJ 70 . Eye phenotypes were scored manually from rank 1 to 10 based on severity, with rank 1 assigned to wild type-like and rank 10 for the most severe phenotype, as described previously 40 . We developed a computational method called Flynotyper (flynotyper.sourceforge.net) that calculates a phenotypic score based on alterations in the hexagonal arrangement of ommatidia in the fly eye 40 . The Flynotyper software detects the center of each ommatidium (orange circle), and calculates the phenotypic score based on the number of ommatidia detected, the lengths of six local vectors with direction pointing from each ommatidium to the neighboring ommatidia, and the angle between these six local vectors ( Figure   3A (i)). Using Flynotyper, we obtained quantitative measures of fly eye roughness with single gene or pairwise gene knockdown. The significance of Flynotyper results compared to a GD control was determined using one-tailed or two-tailed Mann-Whitney tests (Table S9). We found no significant differences in Flynotyper scores between GD and KK control Drosophila lines with and without Dicer2, and therefore we used a single control for statistical analysis ( Figure   S1B). We have previously shown a strong concordance between manual scores and phenotypic scores 40 . In this study, we used manual scoring in conjunction with Flynotyper, as certain features such as necrotic patches, glossy eyes, and overall eye size are not detected by

Confocal Microscopy and image analysis
We acquired Z-stack images of larval eye discs (proliferation assay), pupal eye discs (cellular architecture), and body wall muscles 6 and 7 in the abdominal segments A2 and A3 (NMJ architecture) using an Olympus Fluoview FV1000 laser scanning confocal microscope (Olympus America, Lake Success, NY). Acquisition and processing of images was performed with the Fluoview software (Olympus). We used one or two optical sections for larval and pupal eye disc images, and maximum projections of two or three optical sections were used for NMJ images.
For BrdU staining and proliferation (anti-pH3) assays, maximum projections of all optical sections were generated for display. Area, length, perimeter, and number of branches in neuromuscular synapses were calculated using the Drosophila_NMJ_morphometrics macro in ImageJ 72 . The bouton counts in each NMJ and pH3-positive cells from larval tissues were counted using the Cell Counter Plug-In within ImageJ 70 . We also calculated the number of pH3 positive cells using the Analyze particles function in ImageJ, and found a high correlation (Pearson correlation, r=0.9599, p<0.0001) with counts obtained from Cell Counter Plug-In ( Figure S1F). Significance of cell counts or NMJ features from confocal microscopy compared to GD controls was determined using one-tailed or two-tailed Mann-Whitney tests (Table S9).

Differential expression analysis of transcriptome data
We performed RNA sequencing of samples isolated from fly heads of Elav-GAL4>Dicer2 crosses for MAPK3 rl , KCTD13 CG10465 , DOC2A rph , CORO1A coro , C16ORF53 pa1 , and CDIPT pis , and compared gene expression levels to VDRC control flies carrying the same genetic background. We prepared cDNA libraries for three biological replicates per knockdown model edgeR 76 v.3.16.5 (generalized linear model option) was used to perform differential expression analysis; one biological replicate each for CDIPT pis , CORO1A coro and C16ORF53 pa1 was excluded from the differential expression analysis due to poor clustering with the other replicates. Genes with a log2-fold change greater than 1 or less than -1, and with a corrected false-discovery rate less than 0.05, were considered as differentially expressed (Table S6). We used the log-fold change in expression to confirm reduced gene expression of each 16p11.2 homolog in the tested RNAi lines. These values were similar to expression values obtained by qPCR; we found a positive correlation between qPCR and RNA-Seq derived expression values for 186 differentially expressed genes across the six knockdown models (Pearson correlation, r=0.4258, p=1.372×10 -9 ). Human homologs of differentially-expressed fly genes were identified using DIOPT 60 v.5.3.

Functional enrichment in differentially expressed genes
We used gene set enrichment analysis to summarize the genome-wide list of genes and their levels of differential expression into biological pathways and processes perturbed by knockdown of 16p11.2 homologs ( Figure S7A). First, we used DIOPT 60 to identify fly homologs of all annotated genes in each human Gene Ontology Biological Process term 77 . We then calculated Z scores for all GO terms with less than 500 genes (in order to exclude very general GO terms) across the six knockdown models using the Parametric Analysis of Geneset Enrichment procedure 53 . This method averages the log-fold change in expression of all genes in every GO term, and then subtracts the mean and divides by the standard deviation of the log-fold change levels in all genes. The Z-score represents the degree of up-or down-regulation of all genes within the GO term. We estimated a p-value for each z-score by comparing to the standard normal distribution (two-sided test), and corrected for multiple hypothesis testing using the Benjamini-Hochberg method. 516 GO terms with corrected p-values <0.01 are listed in Table   S6. We also examined enrichment for neurodevelopmental disease [10][11][12][13] and function 78 among differentially-expressed genes ( Figure S7B). We used Cytoscape 79 to visualize the network of cell proliferation (GO:0008283) and cell cycle (GO:0007049) genes that were differentially expressed in the knockdown models ( Figure S7C).

Analysis of 16p11.2 gene interactions in the context of a human brain-specific gene network
We used a human brain-specific gene interaction network 54 to further contextualize the observed interactions in 16p11.2 homologs. This network was built using a Bayesian framework that integrated brain-specific signals from genomic data published in over 14,000 publications, as described previously 54 . Within this network, we mapped 49 genes with identified interactions in the fly eye, and calculated the shortest paths between these genes. This procedure identified 428 additional genes in the network that were critical in connecting the 49 assayed genes to each other (Table S7). We then examined these connector genes for enrichment of genes with cell proliferation and cell cycle GO annotations using a one-sided Fisher's exact test.

DATA AVAILABILITY
Gene expression data for the six 16p11.  We identified 14 homologs of 16p11.2 deletion genes in Drosophila melanogaster (Top), and evaluated global, neurodevelopmental and cellular phenotypes. We also performed transcriptome sequencing and assessed changes in expression of biologically significant genes (Left). Next, we identified modifiers of the one-hit eye phenotype for select homologs using two-hit interaction models. A subset of these interactions was further assessed for cellular phenotypes in the two-hit knockdown eyes. We incorporated all fly interactions into a human brain-specific genetic interaction network (Right).      Figure 2F). (G) A human brain-specific genetic interaction network of all tested 16p11.2 genes and modifier genes, as well as neighboring connector genes. Network nodes with thick borders represent tested genes, with node shape representing gene category. The size of the nodes is proportional to how many connections they have in the network, and the thickness of 36 the edges is proportional to the number of critical paths in the network using that edge. Purple nodes are genes annotated with cell proliferation or cell cycle GO terms.