Identification of a genetic network for an ecologically relevant behavioural phenotype in Drosophila melanogaster

Pupation site choice of Drosophila third‐instar larvae is critical for the survival of individuals, as pupae are exposed to various biotic and abiotic dangers while immobilized during the 3–4 days of metamorphosis. This singular behavioural choice is sensitive to both environmental and genetic factors. Here, we developed a high‐throughput phenotyping approach to assay the variation in pupation height in Drosophila melanogaster, while controlling for possibly confounding factors. We find substantial variation of mean pupation height among sampled natural stocks and we show that the Drosophila Genetic Reference Panel (DGRP) reflects this variation. Using the DGRP stocks for genome‐wide association (GWA) mapping, 16 loci involved in determining pupation height could be resolved. The candidate genes in these loci are enriched for high expression in the larval central nervous system. A genetic network could be constructed from the candidate loci, which places scribble (scrib) at the centre, plus other genes known to be involved in nervous system development, such as Epidermal growth factor receptor (Egfr) and p53. Using gene disruption lines, we could functionally validate several of the initially identified loci, as well as additional loci predicted from network analysis. Our study shows that the combination of high‐throughput phenotyping with a genetic analysis of variation captured from the wild can be used to approach the genetic dissection of an environmentally relevant behavioural phenotype.

choice under laboratory conditions. The suitable location selected for pupation can immediately impact survival in various environmental conditions (Riedl, Riedl, Mackay, & Sokolowski, 2007). For instance, in a humid environment, animals pupating close to the food are less likely to survive, with increased risk of microbial attack or being drowned as other larvae churn up and liquefy the substrate (Markow, 1979;Rodriguez et al., 1992). The pupation site choices in Drosophila melanogaster were revealed to be under stabilizing selection, with significantly decreased viabilities for individuals pupating both close and far from the medium surface (Joshi & Mueller, 1993).

Artificial selection experiments on pupation height in several
Drosophila species revealed a quick response in both directions (i.e., low and high pupation) of selection, with significant divergence being observed after only a few generations of selection (Casares & Carracedo, 1986;Garcia-Flores, Casares, & Carracedo, 1989;Singh & Pandey, 1993). The rapid response of artificial selection resulted in realized heritability estimates ranging from 0.13 to 0.18. Genetic association studies were conducted to explore the genetic basis of pupation site choice in Drosophila and were consistent with it being a complex behavioural trait with a polygenic basis (Bauer & Sokolowski, 1985;Erezyilmaz & Stern, 2013;Riedl et al., 2007;Sokolowski & Bauer, 1989). The genetic contributions were mostly found on the autosomes, but which chromosome contributed most to the explained variation remained disputed (Bauer & Sokolowski, 1985;Sokolowski & Bauer, 1989).
Using a set of 76 recombinant inbred strains for a quantitative trait locus (QTL) study, Riedl et al. (2007) identified a major QTL at position 56A01-C11 on chromosome 2R, which includes at least 39 annotated genes.
Natural populations of D. melanogaster harbour a wealth of genetic variation, which has been proven to constitute a powerful resource to detect the genetic architecture for a wide range of phenotypes (Huang et al., 2014;King, Macdonald, & Long, 2012;MacKay et al., 2012). The Drosophila Genetic Reference Panel (DGRP), a population consisting of more than 200 highly inbred D. melanogaster strains derived from the Raleigh, US population (Huang et al., 2014;MacKay et al., 2012), has been successfully exploited to identify the association between a broad range of phenotypes and their underlying genetic basis, including several behavioural traits (Lee et al., 2017;Rohde, Gaertner, Ward, Sørensen, & Mackay, 2017;Shorter et al., 2015;Xue, Wang, & Zhu, 2017). The highly inbred nature of these DGRP strains allows accurate association mapping, through repeatedly sampling for individuals with the same genotype. Together with the well-resolved richness in genetic polymorphisms and rapid decay in linkage disequilibrium (LD) in these strains (Huang et al., 2014), these make the DGRP an excellent panel to study the genetic basis of pupation height choice in D. melanogaster at a fine-scale resolution.
The efficient and accurate measurement of phenotypes is often a limiting factor in QTL studies, particularly for behavioural traits where variation may be impacted by large numbers of interacting loci. Automated or semi-automated measurements can solve this problem. We have previously described an automated procedure for pupal case size measurements (Reeves & Tautz, 2017). Here, we used the same general approach as a means for the high-throughput accurate measurement of pupation height, while attempting to control for nongenetic confounding factors. We systematically assessed the pupation height choice patterns in both a range of natural global stocks and 198 DGRP inbred strains. We also examined the potential influence of several biotic factors, including sex, Wolbachia infection and parental effects. We explored the genetic heritability of pupation height and dissect the contributions of genetic variants via a GWA analysis. This allowed us to identify 16 significant genetic loci associated with pupation height variation across DGRP strains, of which we experimentally evaluated 7 candidates in gene disruption lines. Furthermore, we also analysed both gene expression data and genetic interaction data, in order to identify a genetic network involved in pupation height choice in D. melanogaster.

| Drosophila strains
The list of wild-type strains, DGRP inbred lines, deficiency and transposon insertion stocks, and their progenitor lines used in this study and the detailed information is provided in Table S4. Flies were reared under standard culture conditions (cornmeal-molasses-agar medium, 24°C, 55%-75% relative humidity, 12-hr light/dark cycle).
A HOBO® data logger was placed in the incubator to monitor and record the environmental changes, that is, temperature, light and humidity, across all the experimental periods.

| Automated phenotyping of pupation height
A previously established automated pupal case length detection pipeline was adopted and modified for the automatic screening of pupation height measurements (Reeves & Tautz, 2017). In brief, standard food was dispensed into 28.5-mm-diameter and 95-mmheight vials (Genesee Scientific), and the food height (defined as the distance from the surface of the food to the bottom of the vial) for each vial was manually measured and recorded. Once the food vials had fully cooled, 10.1 cm × 10.5 cm squares of overhead projector film (Nobo, plain paper copier film, 33638237) were slid into each vial lining their entire vertical wall. Approximately 10 healthy female flies (15 for inbred stocks) and 5 healthy male flies were introduced into each vial, for which a custom-printed semitransparent label (GA International Inc.), including a unique barcode, was affixed to the outside of each vial. Adult flies were removed from the vials after 1-2 days, and vials were kept under the same incubation condition (see above) for another 8-9 days to allow them to reach pupation stage. In general, by the 10th day after the parents were initially introduced, the majority of offspring in the vials were present as pupae attached to the transparent film. The film was gently moved out from each vial, the food from the lower part was scraped away, and any larvae or pupae at white puparium stage (P1) were removed.
The film was then placed into a pre-made plastic frame, which holds the film flat for further photographing using bottom illumination in a light-tight box. Batches of the resulting images were then introduced into the image analysis procedure.
The open-source image analysis software CellProfiler (v2.1.0) (Carpenter et al., 2006) was applied for the simultaneous recognition of pupae and measurements of a variety of attributes (Dataset S1), with a customized pipeline slightly modified from Reeves and Tautz (2017). In brief, any shape significantly distinct from the background was first identified without size restriction (module "identify primary object"). The identified objects composed of multiple juxtaposed pupae were disentangled into distinct pupae (module "untangle worms"). The resulting pupae were then shrunk and re-delineated based on boundary changes in pixel intensity (module "identify secondary objects"). Finally, a specific confidence level was assigned for each pupa based on its size attribute. The digital outlines of pupae were overlaid onto a cropped version of the original image for easy visualization ( Figure 1). A manual check on 40 randomly selected films showed that the CellProfiler pipeline can successfully identify 96% of true pupae (sensitivity), with an accuracy of 81% for identified putative pupae. To further improve the detection accuracy, an additional refinement criterion was defined based on the size attributes of "true" pupae from manual curation. Applying the new criteria, the accuracy for pupae detection was improved to around 99.85%, with only a negligible fraction (<0.7%) loss of true-positive results.
In addition, a 1 euro cent coin (16.25 mm diameter) was included in each image, for the control of camera coordinate changes and the conversion of measurements from pixels to millimetres. The pupation height for each pupa was calculated as the subtraction of the vial food height from the vertical coordinate measurement (CellProfiler parameter: "Areashape_Y"). Overlaid images and files with a variety of attribute measurements were imported into a FileMaker database (v14,FileMaker Inc.). The quality filtering of pupae and related analysis were conducted with the tools implemented in the database. A more detailed description on the pupae detection pipeline can be found in File S1.

| Treatment of confounding factors
Pupal density in the vial is a biotic factor that could affect the pupation site selection preference of third-instar larvae (Joshi & Mueller, 1993;Sokolowski & Hansell, 1983). Here, density was controlled by limiting the number of parents used per vial and restricting the number of nights they remained before being cleared (see above). Further, to reduce the possible bias from including vials containing very low number of individuals, only vials with a pupal density of a minimum of 15 were considered as reliable, and a measurement for each stock should include at least 6 such reliable vial measurements.
All of the tested stocks exhibited a similar positive correlation between individual density and pupation height estimate (File S1). The following equations were used to correct the influence of individual density in the vial on the mean estimate of pupation height in a first step:

| Automatic measurement of pupal case length
The measurement of pupal case length followed the procedure described in Reeves and Tautz (2017). In brief, the pupal case length is defined as the length of the major axis of the ellipse that has the same normalized second central moments as the region of identified pupae, measured with the "Areashape_MajorAxisLength" index in CellProfiler. As the pupal case length measurements are relatively robust to the pupal density in the vial and the minor change of incubator humidity (File S1), the measurement for pupal case length was not corrected for these factors. A reciprocal crossing approach was used to detect whether any maternal, for example, genetic effect from mitochondria, or paternal F I G U R E 1 Automated measurement of pupation height. (a) The larvae in the vial crawl up the wall for pupation. The vial includes a plastic sheet that lines the wall. This is taken out after a sufficient number of pupae are attached, and it is photographed (b). The image analysis software CellProfiler is then used to identify the outlines of the pupae (c) and to identify single ones that can be reliably measured (marked in blue). After optimizing the system, we achieve accuracy for pupae detection of around 99.85%, with only a small fraction of loss of true-positive results (<0.7%) [Colour figure can be viewed at wileyonlinelibrary. com]

| Wolbachia infection, sex and maternal/ paternal effect test
effects for pupation site selection exist. Two pairs of high and low pupation height DGRP inbred lines (DGRP-486 and DGRP-383; DGRP-228 and DGRP_318) were selected for this analysis. The pupation height of F1 offspring from reciprocal crossing, that is, virgin females from high stock crossing with males from low stock, and vice versa, was measured and compared with the phenotype of their parental stocks.

| Estimates of heritability
The broad-sense heritability (H 2 ) was estimated with the variance components of a linear model of the form: Phenotype = Population mean + Line effect + error (Schmidt et al., 2017). Total phenotypic variance was estimated as Genetic Variance + Environmental Variance, and the H 2 was thus estimated as The narrow-sense heritability was estimated as the proportion of variance in a phenotype explained by all available genetic variants used for mapping, an estimate that is often called "SNP heritability" (Visscher, Hill, & Wray, 2008;Wray et al., 2013). In practice, a genetic relationship matrix (GRM) between pairs of inbred strains from all the DGRP annotated genetic variants was built by using GEMMA (version 0.96) (Zhou & Stephens, 2012), and then, the narrow-sense heritability (denoted as PVE) was calculated based on the above GRM with the univariate linear mixed model (Zhou & Stephens, 2012) implemented in GEMMA.

| Principal component analysis (PCA) and GWA analysis
The genetic variant information and major genomic inversion status were retrieved from DGRP freeze 2 (Huang et al., 2014). Genetic variants with missing values above 20% and minor allele frequency below 5% were excluded from further analysis, with 1,903,028 genetic variants (on average ~14 genetic variants per kb) passing the stated criteria. To assess the possible influence of population structure on the pupation site selection, the PCA module from PLINK v1.90 (Purcell et al., 2007) was used to identify top principal components (PCs) from the filtered genetic variant data. The projection length of each strain on top 20 PCs was used to test the influence of cryptic population structure on pupation site selection.
The linear regression model implemented in PLINK was initially used to perform association analysis for the above filtered genetic variants without correction of population structure. In order to correct any potential influence from cryptic population structure, a linear mixed model using the fastlMM (Lippert et al., 2011) program (version 0.2.32) was also applied for GWA analysis. By using fastlMM, another two GWAs were conducted for without and with the incorporation of pupal case length.
A nominal p-value threshold (p-value = 1 × 10 −5 ), commonly used in Drosophila GWA studies (Dembeck et al., 2015;Lee et al., 2017), was applied to identify candidates for genetic variants associating with pupation height (note that this is not meant to be a significance cut-off). The r package "qqMan" (Turner, 2014) was exploited for the visualization of GWA results in a Manhattan plot and QQ plot.
The standardized effect size for each genetic variant was calculated as one-half the difference between marker classes divided by the overall phenotypic standard deviation (Harbison, McCoy, & Mackay, 2013). The associating genes for each genetic variant were annotated by SnpEff (Cingolani et al., 2012) with default parameters. In brief, all the protein-coding genes within 5 kb up/downstream of the target genetic variant were taken as associating genes.
The genotypic linkage disequilibrium (LD) for each pair of significant genetic variants was tested by calculating the squared correlation estimator r 2 with PLINK (Purcell et al., 2007). Moreover, the r 2 values for each genetic variant and all other genetic variants were also computed. A significant genetic region (QTL) was defined by the position of the most distant downstream and upstream genetic variants showing a minimum r 2 of .8 to the significant genetic variants.
All the associating genes as determined above, together with the genes within the LD regions, were considered as candidate genes for further analysis.

| Expression and genetic interaction network analysis
The gene expression profiling data (downloaded in December of 2017) from Drosophila modENCODE project (Brown et al., 2014) were used for expression analysis. These were generated by measuring the genome-wide gene expression for 5 tissues in third-instar larvae stages, including central nervous system (CNS), digestive system, fat body, image disc and saliva glands. The expression level for each gene within each tissue was measured in units of RPKM. The fractions of genes with expression (RPKM > 0), as well as the gene expression levels of both GWA candidate genes and network-predicted genes (as stated below), were compared with those of total annotated protein-coding genes (Fisher's exact test for the comparison of the fractions of expressed genes; Wilcoxon rank-sum test for the comparison of gene expression levels).  (Shannon et al., 2003) was used for the visualization of these subnetworks.

| Functional validation experiments
Functional validation experiments to test the phenotypic effects of  Note that the detected significances of the phenotyping effect of candidate genes from this approach are likely to be underestimated, as in theory only half of the individuals in the experiment group contain the gene semideletion.
Moreover, phenotypic effects of additional four genes, that is, epidermal growth factor receptor (Egfr), E2F transcription factor 1 (E2f1), Ras oncogene at 85D (Ras85D) and p53, in the above-predicted network were further checked, via direct comparisons of pupation height status between the co-isogenic progenitor stock and transposon disruption of each target gene (Table S4D). For gene p53, two different transposon insertion strains (with disruption at different gene regions) were tested. The experimental verification scheme follows the protocol as described in the above text.

| Establishment of an automated phenotyping pipeline
The acquisition of phenotype data from a large number of individuals is a prerequisite for high-resolution genetic mapping studies. Instead of using manual measurement-based approach from previous studies on pupation height (Bauer & Sokolowski, 1985;Erezyilmaz & Stern, 2013;Riedl et al., 2007;Sokolowski & Bauer, 1989), we adapted an image analysis-based phenotyping pipeline (Reeves & Tautz, 2017). This pipeline was initially developed for the high-throughput measurement of pupal case length and was shown to have the capability for the automatic detection of pupae with a high precision (Reeves & Tautz, 2017). We modified it for the purpose of the automated measurement of pupation height, defined as the distance from the vertical coordinate of pupation site (pupal centre) to the food surface in the vial in millimetre (mm). Figure 1 shows an example of the automated measurement of pupation height. The detailed experimental setup and the entire phenotyping procedure are provided in File S1.
Density of individuals within the vial is a common major environmental covariate of many Drosophila traits, including pupation site choice (Joshi & Mueller, 1993;Sokolowski & Hansell, 1983). In the present study, individual density variation was controlled in an indirect manner through limiting the number of parents used per vial (10 females for wild-type strains, and 15 females for inbred strains) and restricting the number of nights they remained before being removed (1-2 nights). Still, some variation was apparent that needed to be addressed. Based on a set of test experiments, we defined a minimal sampling rule of at least 15 measured pupae per vial and at least six replicate vials per strain. We found that almost all stocks exhibited a uniformly positive relationship between individual density and pupation height, allowing us to use an average slope (0.145) across all tested stocks to correct the mean estimate of pupation height in all vials (File S1). However, note that by using such a single correction factor, we ignore the possibility of strain-specific effects, which leads to somewhat increased noise in the data.

| Wild-type variation is captured by the DGRP stocks
Two distinct sets of strains were used to explore the variance in pupation height choice. The first set consisted of 14 natural wild-type D. melanogaster strains, collected from different parts of the world (Table S4A). The second set was from the DGRP (Huang et al., 2014) and included 198 lines (Table S4B). In order to correct for environmental factors, especially cryptic differences in humidity (Casares et al., 1997;Sokal et al., 1960), two wild-type stocks (S-317 and S-314) representing two extremes of pupation height from the first strain set were continually remeasured to act as controls throughout all experiments (File S1). The estimates of pupation height for the strains were corrected based on the average pupation height of the two control stocks across all rounds of experiments. Similar as above, we note that this procedure ignores the possibility of strain-specific effects, that is, is expected to increase the noise. In the following, all  The wild-type lines showed no obvious geographical clustering of pupation heights. The spread of pupation height among the DGRP stocks exceeds that of the wild-type stocks, suggesting that they may capture at least a major part of the existing variation in D. melanogaster.

| Sexual dimorphism and parental effects
Sexual dimorphism, the condition where sexes from the same species exhibit different characteristics for morphological or behavioural traits, is a commonly observed phenomenon (Berner, Sládek, Holt, Niskanen, & Ruff, 2017). Regarding the pupation site choice in Drosophila, a controversy on the existence of sexual dimorphism has persisted for several decades. Early studies have reported no sexual dimorphism (Markow, 1979;Sokolowski & Bauer, 1989;Welbergen & Sokolowski, 1994), while the results from other studies showed that males pupate significantly higher than females (Casares & Carracedo, 1987;Riedl et al., 2007).
To address the sexual dimorphism question, a distinct data set incorporating both pupation height and sex information of >4,000 randomly selected (2,340 females and 1,935 males) individuals from 728 vials generated from the 4-way pedigree data set reported by Reeves and Tautz (2017) was analysed. Since this earlier study had not recorded the level of the food surface, the pupation height for each sexed pupae was calculated as the deviation from the corresponding vial to average pupation height (Figure 3a) (Reeves & Tautz, 2017). As shown in Figure 3b, there was a significant difference in pupation height between males and females (Wilcoxon rank-sum test: p-value = 1.5 × 10 −7 ), with male individuals pupating on average around 2 mm higher than females. This result is roughly consistent with the observed sexual dimorphism reported in two previous studies (Casares & Carracedo, 1987;Riedl et al., 2007).
A parental bias, by which the phenotype of an individual depends more upon the mother's or father's phenotype or genotype, can be observed for some traits (Reik & Walter, 2001). This can be caused by the inheritance of genetic material in the cytoplasm (e.g.,  (Sokolowski & Bauer, 1989).
Wolbachia pipientis is a maternally transmitted endosymbiotic bacterium that infects around 53% of DGRP strains (Huang et al., 2014). It was reported to have a significant effect on some behavioural traits, for example, acute and chronic resistance to oxidative stress (Huang et al., 2014). Two different approaches were used to explore a possible effect of Wolbachia infection on pupation height. First, the statistical analysis on the pupation height between all tested strains with infected and noninfected strains exhibited no significant difference of pupation height (Wilcoxon rank-sum test: p-value = .29, Table S1). Second, the experimental phenotypic comparison between pupation height of three randomly selected DGRP strains with Wolbachia infection and those after the removal of Wolbachia using tetracycline treatment showed no significant statistical difference on pupation height choice for all tested strains ( Figure S1). Accordingly, the Wolbachia infection on DGRP strains was not incorporated in the association analysis below, as both the indirect and direct evidences described above revealed no strong effect on pupation height choice.

| Heritability and chromosome effects of pupation height
The broad-sense heritability (H 2 ) was estimated by determining the proportion of total variance in the mean strain height measurements compared to the average within each strain (Schmidt et al., 2017).
The narrow-sense heritability (h 2 ) was estimated here as "SNP heritability" (Wray et al., 2013), that is, the estimate of the proportion of phenotypic variance explained by all available SNPs (or genetic variants) in the DGRP stocks.
All estimates are shown in Table 1. Values of H 2 of 0.64 (0.70 for wild-type strains) and h 2 of 0.46 (SE: 0.2) based on the estimates from DGRP inbred stocks imply higher heritability than from previous estimates within this species (Casares & Carracedo, 1986;Garcia-Flores et al., 1989;Singh & Pandey, 1993).
Partitioning the variance by chromosome reveals that all chromosomes, except the 4th, contribute a substantial part to the variance of pupation height ( Figure S2). The minimal contribution from chromosome 4 can be ascribed to the limited number of genetic variants within this chromosome. In line with previous reports (Bauer & Sokolowski, 1985;Sokolowski & Bauer, 1989), we find a somewhat higher contribution of autosomes to the variance of pupation height and also a slightly larger effect from chromosome 2 compared to chromosome 3 (Bauer & Sokolowski, 1985). These observations also confirmed the polygenic nature of pupation height choice in D. melanogaster (Bauer & Sokolowski, 1985;Erezyilmaz & Stern, 2013;Riedl et al., 2007;Sokolowski & Bauer, 1989).

| GWA analysis
The GWA analysis was based on the genetic variants of DGRP freeze 2 (Huang et al., 2014); variants with missing values above 20% and minor allele frequency below 5% were excluded from further analysis. A few possible covariates were first assessed, in order to optimize the model used for the GWA analysis.
To investigate whether any cryptic population structure could contribute to the observed variation in pupation site choice of inbred stocks, PLINK (Purcell et al., 2007) was used to identify major PCs of genetic variants in the DGRP strains. Only one major cluster was found in these stocks based on PC analysis (Huang et al., 2014), and there were no obvious clusters of strains that have different pupation height status on the first two principal axes ( Figure S3). Moreover, only three out of the top 20 PCs showed significant (Pearson's correlation test, p-value < .05) correlations with pupation height. Though all these correlation coefficients are very low (ranging from 0.002 to 0.245, Figure S4), the slight correlations seem more likely to come from a few outlier individuals.
Given their driving force on population divergence and speciation (Hoffmann & Rieseberg, 2008), major genomic inversions in DGRP strains (Huang et al., 2014)  two PCs of genomic variation showed significant effects only from In(2L)t and In(3R)Mo (Table S2), indicating their possible roles in population subdivisions (Hoffmann & Rieseberg, 2008). However, there is no significant association between pupation height and all tested genomic inversions, including In(2L)t and In(3R)Mo (Table S2).
Consequently, genomic inversions were not included in the following association mapping analysis.
The original DGRP lines were constructed such that population structure effects should be minimized, but some genetic relatedness leading to cryptic population structure might still exist (Mackay & Huang, 2018). In order to correct any major influence from the cryptic population structure that we identified above, a linear mixed model  TA B L E 1 Statistics for the estimation of the heritability of pupation height choice without correction of population structure, confirming the minimized population structure effects within DGRP lines ( Figure S5). Still, for further analysis, we use the GWA with correction of population structure (Figure 4).
A previous study has shown a possible role of larval size on their pupation site choice (Vandal et al., 2012). Here, we use pupal case length as a proxy of larval size (Reeves & Tautz, 2017). We find that there is indeed a weak, but significant negative correlation between pupal case length and pupation height, though only at vial level ( Figure S6a, Pearson's correlation test, correlation coefficient r value = −.14, p-value = 3.1 × 10 −8 ), but not at strain level When applying a nominal GWA p-value threshold of 1 × 10 −5 (defined below), we found that these two GWAs shared most of the identified significantly associated genetic variants ( Figure S6d).
Based on these analyses, we conclude that the identified candidate genetic variants on pupation height are mostly independent from its association with pupal case length. Therefore, pupal case length was not included in the GWA for further analysis.

| Candidate genes from the GWA analysis
A nominal GWA p-value threshold of 1 × 10 −5 (Dembeck et al., 2015;Lee et al., 2017) was set to detect significant genetic variants associated with pupation height. At this cut-off, we found 16 significant genetic variants (13 SNPs and 3 indels) associated with pupation height in the DGRP strains, corresponding to 38 associating genes that locate within 5 kb up/downstream (default setting in SnpEff (Cingolani et al., 2012)) of these genetic variants (Table 2).
A moderate negative correlation was observed between the minor allele frequencies (MAFs) and their GWA association p-values of these genetic variants (Pearson's correlation test, r value = −.50, p-value = .05), indicating a tendency that gene variants with high MAFs are more likely captured by GWA test (Table 2 and Figure   S7).
To identify additional candidate genes associated with the variants, we examined the long-range linkage disequilibrium (LD) between pairs of detected candidate variants and with other genetic variants found in the DGRP strains. No significant linkage between physically distant (≥1 Mb, with r 2 ≥ .8) genetic variants was found, suggesting the associations are not confounded by long-range LD.
LD blocks were then calculated for each significant genetic variant with a commonly used threshold r 2 = .8 (Pallares, Harr, Turner, & Tautz, 2014), and 6 significant LD blocks were found with average block size of 6.96 kb (Table S3)

| Phenotype confirmations
The advantage of Drosophila as a model system is that one can use mutant alleles that have been constructed in a common coisogenic background to test whether different alleles in genes implicated by the GWA analyses indeed affect pupation site choice.
Despite the potentially low effect of each candidate locus, an observable phenotypic change would still be expected, given the great sample size from our automatic phenotyping system. Gene disruption lines were available for 7 candidate loci (with disruption F I G U R E 4 Manhattan plot for GWA results on pupation height with the correction for population structure. The p-values in −log10 transformation are shown in y-axis. The dashed horizontal line marks the significance p-value (10 −5 ) threshold used in this study. The genetic variants in rectangles are the ones applied for experimental verification in Figure 5 in the gene region of at least 1 gene for each locus) within the 16 identified associated loci. Two constitute transposon insertion mutations (Bellen et al., 2011), and five constitute small deficiencies (Table S4C). Experimental tests involve replicated phenotypic comparisons between the co-isogenic progenitor stock and homozygous or heterozygous disruptions of the target genes (see Methods).
An overview of the measured pupation height differences compared to the respective progenitor stocks is provided in Figure 5.
Five out of 7 tested candidate loci showed a significant difference, two with an increased height and three with a decreased one. Both of the two transposon insertion stocks that are homozygous viable (Smrter -Smr, scrib) show a decreased pupation height. One of these two insertion stocks (scrib) reached the statistical significance (pvalue < .05), while the other one (Smr) did not (p-value = .21).
The other stock that shows no significant overall change (Arl6IP1/nonA-l/Fas1) is a deletion stock that is homozygous lethal, that is, the fact that we did not find an effect on pupation height

| Expression and network analysis
Pupation follows the third-instar larval (L3) stage in Drosophila.
Hence, one can ask whether genes involved in pupation site choice are more likely to have a substantial expression at this stage, especially in the central nervous system (CNS), given that this is a behavioural phenotype and such a CNS enrichment is also known from another behavioural phenotype (Ayroles et al., 2015). To test this hypothesis, the gene expression profiling data from the Drosophila modENCODE project (Brown et al., 2014) were explored, which measure the genome-wide gene expression for five tissues in L3 stages in D. melanogaster, including CNS, digestive system, fat body, imaginal disc and salivary glands. Larger fractions of expressed genes (RPKM > 0) for the identified candidate gene data set were

TA B L E 2 GWA results on pupation height
found for most of the tested tissues in L3 stages (including CNS, digestive system and fat body), compared with the expression profiling of total annotated protein-coding genes as controls, but only at moderate p-value levels ( Figure S9a). We also find that expression levels of candidate genes compared to all the annotated proteincoding genes are significantly higher in the CNS (Wilcoxon rank-sum test, p-value = .05), but not other tissues ( Figure S9b). This observation is consistent with the above prediction that the identified candidate genes for pupation site choice are enriched in the CNS of third-instar larvae, where they could have a direct influence on behaviour.
For 13 out of 49 (27%) genes in the 16 candidate loci, their genetic interactions have been documented in Flybase v6.19 (http://www. flyba se.org) (Attrill et al., 2016). We used this information to construct a computationally predicted network of genetically interacting genes, allowing one intermediate gene (i.e., a noncandidate gene connecting two candidate genes). This analysis revealed a connected subnetwork of 5 candidate genes from the GWA analyses and 4 computationally recruited intermediate genes (Figure 6a). The probability that this connected subnetwork would have arisen when the same number of genes is randomly sampled is very low (p < 2.2 × 10 −16 ). Based on the availability of gene disruption stocks, the phenotypic influence of one intermediate gene (i.e., Ras85D) in the above network was experimentally tested and showed a significant effect on pupation height (Table 3). Therefore, this network appears to represent a core network controlling the pupation site choice in D. melanogaster.
To further substantiate this inference, we predicted that the genes that are directly connected to GWA candidate genes may also be involved in the pupation height choice phenotype. Hence, we used the genetic interaction information to construct another computationally predicted network (Figure 6b), including both the above core subnetwork and the first-degree interacting genes of identified candidate genes within the core subnetwork. The phenotypic effects of three randomly selected first-degree interacting genes (Egfr, E2F transcription factor 1 -E2f1 and p53) in the above network were then analysed, via direct comparisons of pupation height status between the co-isogenic progenitor stock and the corresponding gene disruption lines. All these three genes showed indeed a significant phenotypic effect on pupation height (p-value < .05), and the results were consistent for the disruption of different alleles for the same gene (i.e., p53, see Table 3). Gene ontology enrichment analysis of the genes in the network showed significant enrichment for genes associated with cellular and developmental processes, including larvae development and neuron development (Table S5). Therefore, the above network may represent a core network controlling the pupation height choice in D. melanogaster via the development of neuronal connections.

| D ISCUSS I ON
Pupation site choice is an ecology-related behavioural trait that is critical to the fitness of holometabolous insects, including Drosophila (Joshi & Mueller, 1993;Rodriguez et al., 1992;Sokolowski, 1985). Previous artificial selection experiments on pupation height in Drosophila revealed that this trait is heritable (Garcia-Flores et al., 1989;Markow, 1979). However, the genetic basis of pupation height remained largely unexplored, mainly due to the lack of large numbers of reliable measurements. To tackle this issue, we developed a semi-automated phenotyping pipeline to reduce the variance of pupation height measurements, via repeated measurements for a given Drosophila strain in a highthroughput manner.
Applying the semi-automated phenotyping pipeline, we found a substantial variation on pupation height for both a collection of natural stocks and DGRP inbred lines (Figure 2). It is noteworthy that the natural stocks show a smaller range than the DGRP stocks, F I G U R E 5 Phenotypic effects of candidate genes from gene disruption tests. The phenotypic effect was measurement as the deviation of pupation height of stocks with gene disruption compared with that from progenitor stocks via Wilcoxon rank-sum test. The error bars show the standard error of mean values. **: p-value < .01; *: p-value < .05; Stocks marked with "#" were small deficiency stock tested as heterozygotes; the remaining two stocks are transposon insertion stocks tested as homozygotes The existence of sexual dimorphism on pupation height in Drosophila was confirmed, with the observation that male individuals pupate on average ~2 mm higher than females (Figure 3). The difference in pupation site choice between sexes may be due to their distinct developmental timing, as females generally pupate later than males, and later larvae tend to pupate lower, possibly due to a response to diminishing levels of humidity inside the vials (Casares & Carracedo, 1987 (Figure 3), suggesting an additive model of inheritance on pupation site selection in D. melanogaster (Sokolowski & Bauer, 1989).
Our present study allows a more refined systematic exploration of heritability of pupation height in D. melanogaster using alternative methods to the existing artificial selection-based ones F I G U R E 6 Gene interaction subnetwork connecting identified candidate genes. (a) Gene interaction subnetwork that connects 5 candidate genes (circles) through other 4 intermediate genes (rectangles); (b) gene interaction subnetwork that includes both the above subnetwork and the first-degree interacting genes of candidate genes. Highlighted in red are the ones for which the phenotypic effects were experimentally tested (Garcia-Flores et al., 1989;Markow, 1979). Values with H 2 of 0.64 (0.70 for natural stocks) and h 2 of 0.46 (SE: 0.2) based on the estimates from DGRP inbred stocks (Table 1) imply higher heritability than that (realized heritability h 2 = 0.13-0.18) from previous estimates within this species (Garcia-Flores et al., 1989;Markow, 1979 (Shorter et al., 2015), DDT resistance (H 2 = 0.8) (Schmidt et al., 2017) and adult foraging behaviour (h 2 = 0.52) (Lee et al., 2017).
The substantial phenotypic variation on pupation height, as well as the richness of genetic variation in DGRP inbred lines, provides the possibility to dissect the genetic architecture of this behaviour trait at a fine scale. GWA analysis returns 16 significant associating genetic loci (corresponding to 49 candidate genes) from all major chromosomes (Table 2). Expression analysis of these candidate genes shows their enrichment in CNS of third-instar larvae, suggesting the role of controlling pupation height choice through the central nervous system ( Figure S9).
The phenotypic effects on pupation height from several candidate genes were experimental confirmed with gene disruption assays ( Figure 5). These include a potential cation transport function (CG15270), transmembrane transporter activity (Oatp74D), a scaffolding protein regulating apicobasal polarity (scrib) and genes of as yet unclear function (CG7567 and CG7029). scrib was found to be widespread expressed in olfactory organs and the central nervous system (Ganguly, Mackay, & Anholt, 2003), and has been reported to be associated with several behavioural traits in D. melanogaster, including olfactory behaviour (Ganguly et al., 2003), adult foraging (Lee et al., 2017) and sleep (Harbison et al., 2013).
More intriguingly, we predicted a genetic interaction subnetwork based on the connections of candidate genes, with the well-studied gene scrib as a hub gene ( Figure 6). The phenotypic effects on pupation height of four randomly selected genes from the predicted network were further experimentally confirmed (Table 3). These include another well-studied gene Egfr, which is the transmembrane tyrosine kinase receptor for signalling ligands in the TGFα family and was also found to function in neuronal development and behavioural traits in Drosophila (King, Eddison, Kaun, & Heberlein, 2014;Potdar & Sheeba, 2013). Ras85D encodes a protein that acts downstream of several cell signals, most notably of receptor tyrosine kinases (RTKs), and has been reported to be involved in pupal size determination (Li, Li, Jin, Chen, & Ma, 2017). Another component of the network is p53, which is a general regulator of the cell cycle, but which has also been found to be involved in central nervous system development in Drosophila (Bauer et al., 2007) and behavioural traits, such as the entrainment TA B L E 3 Phenotypic effects of core network genes The gene expression profiling data were from Drosophila modENCODE project (Brown et al., 2014 This gene is also a significant associated gene from GWA results. of the circadian rhythm in mice (Hamada, Niki, & Ishida, 2014).
Furthermore, both gene expression analysis and GO functional enrichment analysis on the total gene set from the subnetwork revealed their significant association with third-instar larvae development/pupal morphogenesis and neuro-related functions (Table   S5). Thus, we propose this scrib-centred subnetwork involved in the development of the nervous system to be at least part of a core network that also govern pupation height choice in D. melanogaster. Interestingly, although scrib is inferred as central gene, the actual phenotypic effect of the gene disruption line is not very high ( Figure 5). This suggests that the transposon insertion of this stock has only a weak effect on the phenotype. Note that amorphic alleles of scrib are pupal lethal.

| CON CLUS IONS
In contrast to a previous QTL study on pupation site choice (Riedl et al., 2007), we do not simply find one major locus, but show a polygenic architecture that underlies variation in the pupation height in D. melanogaster. All of the major chromosome arms contribute to the trait, suggesting a broad distribution of influencing genes across the genome. Using a nominal GWA threshold p-value lower than 10 -5 , we find 16 loci (corresponding to 49 candidate genes), of which >70% could be verified by gene disruption experiments. Most of these candidate genes are enriched in CNS of third-instar larvae, indicating their involvement in nervous system development. Our results predict a scrib-centred subnetwork as a core network that might potentially govern pupation height choice in D. melanogaster through controlling neuronal connections during larval development.

ACK N OWLED G EM ENTS
We are grateful to Wen Huang for his help on Drosophila Genetic Reference Panel (DGRP) stocks. We thank the laboratory members for helpful discussions and suggestions. We thank Anita Moeller and Elke Blohm Sievers for their excellent technical help and advice with conducting this experiment. This work was supported by institutional funding through the Max Planck Society.

AUTH O R CO NTR I B UTI O N S
D.T., G.R. and W.Z. jointly conceived and planned the experiments.
W.Z. carried out the experiments and analysed the results, with the support from G.R. All these three authors contributed to the final version of the manuscript. D.T. supervised the project.