Nonlinear expression patterns and multiple shifts in gene network interactions underlie robust phenotypic change in Drosophila melanogaster selected for night sleep duration

All but the simplest phenotypes are believed to result from interactions between two or more genes forming complex networks of gene regulation. Sleep is a complex trait known to depend on the system of feedback loops of the circadian clock, and on many other genes; however, the main components regulating the phenotype and how they interact remain an unsolved puzzle. Genomic and transcriptomic data may well provide part of the answer, but a full account requires a suitable quantitative framework. Here we conducted an artificial selection experiment for sleep duration with RNA-seq data acquired each generation. The phenotypic results are robust across replicates and previous experiments, and the transcription data provides a high-resolution, time-course data set for the evolution of sleep-related gene expression. In addition to a Hierarchical Generalized Linear Model analysis of differential expression that accounts for experimental replicates we develop a flexible Gaussian Process model that estimates interactions between genes. 145 gene pairs are found to have interactions that are different from controls. Our method appears to be not only more specific than standard correlation metrics but also more sensitive, finding correlations not significant by other methods. Statistical predictions were compared to experimental data from public databases on gene interactions. Mutations of candidate genes implicated by our results affected night sleep, and gene expression profiles largely met predicted gene-gene interactions.Author summary: Understanding the molecular bases of phenotypes remains a challenge of complex trait biology. We used a combination of selective breeding, RNA-Seq, and Gaussian Process modeling to determine whether de novo gene expression networks could be derived for sleep duration in Drosophila. We bred flies with long and short sleep times, and sequenced RNA from the flies at each generation of selection. Using a hierarchical Bayesian Generalized Linear Model, we identified genes with altered expression across generation in the selected populations. Gene expression trajectories were largely non-linear across time, however, so we developed a Gaussian Process method to more accurately model the data. The Gaussian Process provides an adaptable framework that adjusts to the complexity of gene expression patterns we observed, eliminating the need to specify or assume a specific polynomial model. The Gaussian Process also enabled us to compute covariances among pairs of genes, elucidating gene expression networks for sleep duration. Follow-up mutational analyses validated the candidate genes’ effects on sleep duration and transcriptional analyses of the mutations largely confirmed gene expression network predictions. The Gaussian Process framework is broadly applicable to gene expression data collected across time.