Abstract

Knowledge of the genetic basis underlying variation in response to environmental exposures or treatments is important in many research areas. For example, knowing the set of causal genetic variants for drug responses could revolutionize personalized medicine. We used Drosophila melanogaster to investigate the genetic signature underlying behavioral variability in response to methylphenidate (MPH), a drug used in the treatment of attention-deficit/hyperactivity disorder. We exposed a wild-type D. melanogaster population to MPH and a control treatment, and observed an increase in locomotor activity in MPH-exposed individuals. Whole-genome transcriptomic analyses revealed that the behavioral response to MPH was associated with abundant gene expression alterations. To confirm these patterns in a different genetic background and to further advance knowledge on the genetic signature of drug response variability, we used a system of inbred lines, the Drosophila Genetic Reference Panel (DGRP). Based on the DGRP, we showed that the behavioral response to MPH was strongly genotype-dependent. Using an integrative genomic approach, we incorporated known gene interactions into the genomic analyses of the DGRP, and identified putative candidate genes for variability in drug response. We successfully validated 71% of the investigated candidate genes by gene expression knockdown. Furthermore, we showed that MPH has cross-generational behavioral and transcriptomic effects. Our findings establish a foundation for understanding the genetic mechanisms driving genotype-specific responses to medical treatment, and highlight the opportunities that integrative genomic approaches have in optimizing medical treatment of complex diseases.

UNDERSTANDING individual variability in response to treatment is one of the fundamental challenges in biological sciences. Treatment can be defined as any expected or unanticipated biotic or abiotic exposure; that is, any change in the environment of the organism, whether in the surrounding environment, such as temperature changes, or in the intrinsic environment, such as metabolic alterations caused by ingestion of chemicals. For many scientific fields—including medicine, animal and plant breeding, and evolutionary biology—knowledge of the mechanisms that underlie variability in how an individual responds to such treatment is important. This could, for example, involve predicting the consequences of climate change on species distribution, obtaining consistent output in livestock and plant breeding across heterogeneous production systems, or providing optimal treatment of conditions that require medical attention. The latter example was specifically investigated in the present study.

The majority of studies that have investigated the genetic contributions to the risk of developing ADHD have focused on dopaminergic and serotonergic systems (Faraone et al. 2005; Li et al. 2006; Gizer et al. 2009). Recently, 12 independent risk loci for ADHD were identified (Demontis et al. 2019), some of which are linked to the regulation of neurotransmitter homeostasis by affecting dopamine levels in synapses. The regulation of dopamine levels seems important since many drugs used to treat ADHD have dopaminergic targets, such as MPH and amphetamine, which increase the concentration of dopamine by blocking dopamine transporters (Volkow et al. 2001). Compared to a known ADHD mouse model, the SorcS2 knockdown line, where decreased activity was observed upon exposure to amphetamine, the activity level in a wild-type (WT) mouse was increased (Glerup et al. 2014). In Drosophila, it has been shown in a memory mutant that the key phenotypes of this mutant, namely attention deficit and hyperactivity, could be rescued with MPH (van Swinderen and Brembs 2010) and, more recently, that night hyperactivity and sleep deficit could also be rescued with MPH (van der Voet et al. 2015).

D. melanogaster is an excellent model system, although we acknowledge some limitations for studying human neurodevelopmental disorders using this species (Pandey and Nichols 2011; van Alphen and van Swinderen 2013; Ugur et al. 2016). The use of D. melanogaster as a model species to study human mental disorders is supported by the fact that important genes implicated in many human diseases are remarkably evolutionarily conserved across species (Pandey and Nichols 2011; Ugur et al. 2016). The central nervous systems of mammals and D. melanogaster are derived from a common evolutionary origin (Hirth and Reichert 1999). This is reflected in many structural and biochemical similarities, e.g., synapses between neurons have common protein architectures and the neurotransmitter compounds—such as acetylcholine, glutamate, γ-aminobutyric acid, dopamine, and serotonin—are conserved in mammals and fruit flies. Many of the characteristics of ADHD, such as impulsivity, inattention, or hyperactivity, can be observed in D. melanogaster. Using D. melanogaster and behavioral traits to specifically study ADHD is supported by studies showing that knockdown of D. melanogaster homologs of candidate genes related to ADHD in humans result in changes in behavior and locomotor activity, and that these changes can be ameliorated by treatment with ADHD medications (van Swinderen and Brembs 2010; van der Voet et al. 2015; Rohde et al. 2016b). In summary, there is solid support for D. melanogaster being a useful model organism to study neurological disorders, including ADHD, and that D. melanogaster exhibits behavioral traits, like locomotor activity, that are relevant for studying ADHD (van Swinderen and Brembs 2010; van Alphen and van Swinderen 2013; Wangler et al. 2015).

To obtain knowledge of the genetic signature underlying variability in drug response, we exposed a WT mass population of D. melanogaster to either MPH or a control treatment to assess the general effect of MPH on a common behavioral trait, namely locomotor activity. We then assessed the effect of MPH treatment on gene expression levels in the same WT population to identify transcripts directly affected by MPH. To confirm the patterns obtained using the WT population, we made use of a system of inbred lines, the Drosophila Genetic Reference Panel (DGRP) (Mackay et al. 2012; Huang et al. 2014). The DGRP allowed us to investigate: (1) whether the set of differentially expressed transcripts identified from the WT mass population was also important for the behavioral response in another genetic background (i.e., the DGRP) and (2) whether the behavioral response to MPH treatment was genetically variable. The DGRP lines have been genotyped for millions of single-nucleotide polymorphisms (SNPs) (Mackay et al. 2012; Huang et al. 2014), which enables the identification of sets of genetic markers that influence trait variability. The findings obtained from the genomic analyses of the DGRP were validated by knockdown of gene expression of putative candidate genes, and furthermore by following the behavioral effects and gene expression alterations induced by MPH across one generation of breeding in the WT mass population.

Materials and Methods

Drosophila populations and breeding

A laboratory WT mass population of D. melanogaster was established from a natural population collected in Odder (55°56′42.46″N, 10°12′45.31″E), Denmark in 2016. The population was established from the offspring of ∼200 inseminated females. The inseminated females produced offspring in individual vials, and five virgin males and five virgin females from each inseminated female were used to establish the mass-bred population. After establishment, the mass-bred population was maintained at 23° under a 12 hr:12 hr light:dark photoperiod with a population size of ∼2500 individuals for two generations before the experiments were initiated.

Phenotypic assay

The phenotypic assay consisted of two parts. First, the flies were exposed to one of the two treatments [control or sucrose (SUC)] or MPH using the Capillary Feeder (CAFE) assay (Ja et al. 2007) for ∼24 hr, after which the flies were relocated from the CAFE assay to the activity chambers. Both assays are described in detail below. In Figure S1 in File S1 an overview of all the experimental parts is presented. In all experiments, the flies being tested constituted a random sample collected across a number of rearing vials (for details see sections in Materials and Methods on “WT males,” “DGRP,” and “RNAi knockdown lines”). We did not keep track of which flies were collected from the individual rearing vials and therefore we cannot account for vial (environment) effects in our statistical models. We acknowledge that this is a limitation of our study but given that a large number of flies (∼12.000 individuals) collected across a large number of vials (∼1.000 vials) were tested, we argue that environmental variation was distributed randomly and that nonindependence of data is therefore not a serious issue in our design.

Capillary feeder:

Flies of age 12 ± 12 hr were sorted by sex under light CO2 sedation and males were transferred to new food vials to recover from the CO2 anesthesia (40 ± 2 hr) before exposure to one of the two treatments. The control treatment consisted of a 5% SUC solution (with 5% green food dye) and the MPH treatment was a 5% SUC solution (with 5% green food dye) containing 1.5 mg ml−1 MPH (Sigma [Sigma Chemical], St. Louis, MO) (see Figure S2 in File S1 for a dose response curve). Commonly, psychostimulants are fed to adult D. melanogaster by supplementing the drug into regular Drosophila medium (Andretic et al. 2005; van Swinderen and Flores 2006; van Swinderen and Brembs 2010; van der Voet et al. 2015); however, a visual control for actual ingestion is lacking using such a feeding procedure, thus we used the CAFE assay for visual control of food intake. Up to seven (2–4-days old) males were transferred to the feeding vials containing two 5-μl capillary tubes extending down into the vial. To minimize evaporation from the capillary tubes, the feeding vials were contained within a tightly sealed container with high humidity (> 90%). Flies were allowed to feed from the capillary tubes for ∼24 hr, after which they were moved to the activity chambers.

Activity assay:

Activity was quantified as the distance moved during a 10-min trial in a circular arena as previously described (Rohde et al. 2016b, 2018). Each activity plate contained 36 circular arenas (six-by-six; diameter: 16 mm; and height: 6 mm), which was placed on top of a light box and enclosed within a separate box to minimize external disturbance. Using an iPad Air 2 (Apple, Cupertino, CA), videos were recorded during a 2–5-hr period after light was turned on (from 09:00 am to 12:00 noon), and the total distance covered during the 10-min trial was computed with the tracking software EthoVision XT (v.10) from Noldus (Wageningen, The Netherlands).

Obtaining and processing locomotor activity data

Data on locomotor activity were obtained, processed, and analyzed for different D. melanogaster populations/lines. We obtained activity data for a WT mass-bred population (Figure S1A in File S1), a large number of DGRP lines (Figure S1B in File S1), and for RNAi knockdown lines. We also assessed cross-generational phenotypic effects of MPH using the WT mass-bred population (Figure S1C in File S1). In the following sections, each of these experiments is described in detail.

WT males:

Quantification of the behavioral effect of MPH on WT males was performed using 72 males exposed to SUC and 72 males exposed to MPH (Figure S1A in File S1). For each treatment type, we controlled the fly density during development by having 20 eggs in each of 20 vials with 7 ml standard medium. Flies used for the behavioral experiments constituted a random sample from the 20 vials per treatment to account for potential vial effects.

The flies were distributed to four activity plates (18 flies per treatment type per plate), and activity was quantified as total distance covered during the 10-min trial. Because both treatments were represented on all activity plates, we adjusted the data for experimental plate effects by correcting all observations by the control treatment on one plate, such that the mean activity within an activity plate of the control flies was equal across all plates: , where is the raw locomotor activity measure for any fly, is the mean of the controls on the plate that all observations were standardized according to (we randomly chose the first plate), and is the mean of the control lines on the activity plate that the observations were from. Observations with a median absolute deviation > 2 within a treatment were defined as outliers and disregarded (Leys et al. 2013) (3/144 observations were regarded as outliers).

The statistical effect of treatment was determined by ANOVA: , where y is the vector of adjusted phenotypic observations, T is a vector containing two treatment groups, and e is the residual term. Statistical significance was obtained as an F-test between the full model and a model neglecting the treatment term .

DGRP:

In total, we assayed 172 DGRP lines with ∼30 males per DGRP line per treatment for the behavioral assay (File S1B in File S1). Density during development of experimental flies was controlled by allowing ∼10 flies to reproduce in vials with 7 ml standard medium for 12 hr, after which they were transferred to a new set of vials with 7 ml standard medium. This was repeated four times in total. Male flies used for testing were collected randomly from the available vials to account for potential vial-specific effects.

The behavioral phenotypes were obtained in blocks over 46 days with 10–400 individual flies tested per day (1–12 activity plates). The number of DGRP lines represented on each activity plate varied, but in all cases both treatments were represented for the lines tested on a particular activity plate.

The behavioral phenotypic values were adjusted for known segregating chromosomal inversions, Wolbachia infection status (information available at http://dgrp2.gnets.ncsu.edu, only 2L_t and 3R_p were significantly associated, Table S2 in File S1), for plate and experimental test day effects, and a positional effect (an indicator variable) to account for observations that were obtained from arenas positioned at the border of the activity plate (Table S3 in FIle S1, this was only done for the DGRP experiments where we observed an effect). For the DGRP lines, we also quantified the average amount of SUC and MPH solution ingested, to investigate if the variation in amount ingested had an effect on the behavioral phenotype. We found no overall effect of the amount ingested (Figure S3 and Table S3 in File S1); therefore, we did not assess consumption for the other experiments. The adjustment was obtained by fitting a linear mixed model including all fixed effects;(1)where y is a vector of repeated phenotypic observations, X and Z are design matrices linking fixed and random effects to the phenotype, b is a vector of the fixed effects, g is a vector of random genetic effects that were assumed to be normally distributed with mean of 0 and variance given by , and e is a vector of residuals assumed to be normally distributed with mean of 0 and variance given by , where I is an identity matrix. The additive genomic relationship matrix, G , was computed using all segregating polymorphic markers [i.e., 1,725,755 at minor allele frequency (MAF) > 0.05] as (VanRaden 2008), where m is the total number of SNPs, and W is a centered and scaled genotype matrix, where each column vector is , is the allele frequency of the ith SNP, and is the ith column vector of the allele count matrix, A , which contains the genotypes coded as 0 or 2 counting the number of the minor allele (genotypes are available at http://dgrp2.gnets.ncsu.edu).

The estimated genetic effects has the dimension of the number of DGRP lines, whereas the residuals has the dimension of the number of observations. Thus, the adjusted phenotype for individual i was obtained as,(2)such that the dimension of matches that of y. The adjustment of phenotypes was done within each treatment.

RNAi knockdown lines:

To reduce gene expression ubiquitously for each of 36 selected genes (identified in the genomic analysis of DGRP lines), we crossed virgin females (collected a maximum of 8 hr after eclosion) from the UAS-RNAi lines with males from the tubulin-GAL4 line (y1w*, P{tubP-Gal4}LL7/TM3, Sb1) (O’Donnell et al. 1994). In addition, we crossed males from the tubulin-GAL4 line with the virgin females (collected a maximum of 8 hr after eclosion) from the UAS-RNAi host lines (Table S1 in File S1) to serve as control lines with normal gene expression. The flies used in these crosses were all density controlled during development (∼20 eggs per vial with 7 ml standard food). For each cross we set up five vials (with 7 ml standard food), each with 10 males and 10 females. After 12 hr, the flies were transferred to new vials and this was repeated four times in total; thus generating 20 vials per knockdown or control line.

Of the 36 UAS-GAL4 crosses, 14 of them produced viable offspring (Table S1 in File S1) that were used in the behavioral assay. To account for potential environmental differences between vials within crosses, we randomly collected flies from the 20 vials. Thus, the flies used in the assays consisted of a random sample from the entire population of individuals from each cross.

The activity data were obtained by having the control line and two UAS-GAL4 lines (both treatments) represented on each activity plate; thus, up to six individuals per genotype and treatment. The average number of individuals tested per UAS-GAL4 genotype per treatment was 26. Because the corresponding control line (i.e., offspring from the host line crossed to the tub-GAL4 line) was represented on all activity plates (ntotal = 226/treatment), we adjusted all observations by the control line using the same procedure as described for the WT males. Observations with a median absolute deviation > 2 within line and treatment were defined as outliers, and disregarded (Leys et al. 2013) (162/1406 observations were regarded as outliers).

To assess statistical differences in responses to treatment between the control line and a given RNAi knockdown line, a two-way factorial ANOVA was fitted:, where y is a vector of adjusted phenotypic observations, L is an indicator vector of the control and UAS-GAL4 lines, T is the treatment effect, is an interaction term between genotype and treatment, and e is the residual term. Statistical significance was obtained as an F-test between the full model and a model neglecting the interaction term.

WT F1 males:

One of the findings from the genomic analyses of the DGRP was that a large fraction of the suggested candidate genes (see description further down) were involved in histone-modifying processes. Hence, this also points toward potential cross-generational behavioral effects of MPH (Vassoler and Sadri-Vakii 2014). To test if MPH led to cross-generational effects, we set up several crosses from the WT mass population (File S1C in File S1). Virgin males and virgin females were collected from the WT mass population (a maximum of 8 hr after eclosion). These flies were grouped into four sets, with each set constituting seven virgin males and seven virgin females. Each set of males and females were exposed to one of four treatment types for ∼24 hr: both sexes on control treatment (SUC), both sexes on MPH treatment, females on MPH and males on control treatment, or males on MPH and females on control treatment. Hereafter, males and females (matched by cross type, i.e., ♀SUC ×♂SUC, ♀MPH ×♂MPH, ♀SUC ×♂MPH, or ♀MPH ×♂SUC) were transferred to new vials containing 7 ml standard Drosophila medium. For three (consecutive) days, flies were transferred to new food vials (i.e., three sets per cross type).

For the activity assay, ∼72 untreated 2–4-day-old F1 males per cross type were used (hereafter designated WT F1). To account for potential differences between the vials, we randomly selected F1 males from the three sets of vials (within a cross type); thus, the individuals used for behavioral assessment constituted a random sample of the population of F1 males from each cross type.

The F1 males (untreated) were distributed on eight activity plates (nine flies per cross type), and activity was then quantified as total distance covered during the 10-min trial. Because offspring from all cross types were present on all activity plates, we adjusted the data for experimental plate effect by correcting all observations by the offspring from the control cross on one plate using the same approach as described for the WT males. Observations with a median absolute deviation > 2 within cross type were defined as outliers and disregarded (Leys et al. 2013) (26/281 observations were regarded as outliers).

Cross-generational effects of MPH on locomotor activity were determined by comparing activity levels of the WT F1 males from the ♀SUC ×♂SUC cross to F1 males from the three other cross types, where one or both parents were exposed to MPH (♀MPH×♂MPH, ♀SUC×♂MPH, or ♀MPH×♂SUC). Whether cross-generational effects were observed was determined by ANOVA: , where y is the vector of adjusted phenotypic observations, T is a vector containing the two contrasting cross types, and e is the residual term. Statistical significance was obtained as an F-test between the full model and a model neglecting the treatment term .

Whole-genome transcription profiling

Immediately after obtaining the measurements of locomotor activities of males from the WT mass population (Figure S1A in File S1) and WT F1 males (♀SUC ×♂SUC and ♀MPH×♂SUC, Figure S1C in File S1), the flies were flash-frozen in liquid nitrogen and stored in a −80° freezer for microarray analysis. Gene expression profiling [GeneChip Affymetrix (Santa Clara, CA) Drosophila Genome 2.0 Array, performed by Eurofins Genomics A/S, Galten, Denmark] was obtained using WT males from both treatments and the WT F1 males from the cross between individuals from the control treatment (♀SUC ×♂SUC), and from the cross between females from the MPH treatment and males from the control treatment (♀MPH×♂SUC) (these individuals were chosen because here we observed the numerically largest difference in activity, Figure 3A).

Transcription profiling was performed in three replicates of 8 flies randomly sampled from the ∼72 males per cross type. The GeneChip Affymetrix Drosophila Genome 2.0 Array contains 18,800 probes, which measure the expression of ∼18,500 transcripts. RNA extraction and preparation of cDNA were performed as described by Dyrskjøt et al. (2003). First, 15 µg of cDNA was fragmented and loaded onto the Affymetrix probe array cartridge. Incubation, washing, and staining were performed in the Affymetrix Fluidics Station 450. The probe arrays were scanned at 560 nm using a confocal laser-scanning microscope (Affymetrix Scanner 3000 7G).

Prior to the analysis of the gene expression array, we performed quality control of the array data using the simpleaffy package for R (Miller 2018). Based on the “Data Analysis Fundamentals” from Affymetrix, sample quality was judged based on average background values, the percentages of genes present, and the ratio of the internal control (actin 3′/actin 5′). All samples were equal in their quality values, and therefore kept for further analyses. The transcriptomic data were normalized using the Robust Multi-array Average (RMA) algorithm (Irizarry et al. 2003) as implemented in the Affy package for R (Gautier et al. 2004). All transcripts were annotated using the Affymetrix annotation file (genome version 3, build 36).

Differentially expressed transcripts were identified using the R package limma (Ritchie et al. 2015) by fitting linear models and empirical Bayesian methods, following a paired Student’s t-test and adjustment for multiple testing by the Benjamini–Hochberg (Padj) approach (Benjamini and Hochberg 1995). The significance level of differential expression of transcripts for the WT males was set to Padj < 0.001, and for the WT F1 males we used a less-stringent cutoff of Padj < 0.05. Because of different significance thresholds in the two experiments, we used another approach to compare the expression profiles of the WT males and the WT F1 males. We used B-statistic, which is the log odds that a transcript is differentially expressed. The probability that a transcript is differentially expressed can be computed as (Ritchie et al. 2015). Here, we used as a cutoff.

Quantitative genomic analyses

To investigate the genetic basis underlying treatment response, we applied an integrative quantitative genomic approach using the genomic feature models (GFM) (Edwards et al. 2016; Rohde et al. 2016a, 2017, 2018, 2019; Sarup et al. 2016; Fang et al. 2017). The GFM approach combines genetic data with other types of biological data, for example known pathways or molecular phenotypes such as gene expression data, to infer the contribution of SNPs located within sets of SNPs defined using the external data—genomic features—on the phenotype. The main advantage of the GFM approach compared to, e.g., traditional genome-wide association studies, besides the ability to utilize and leverage biological information obtained from different sources, is that this method considers the joint contribution from multiple SNPs, which typically would have effect sizes that are too small to be classified as associated markers, unless the sample size is very large.

In the following sections, the methods are described in greater detail. The first section describes the generation of genomic feature sets, the second section contains details on the gene set analysis applied to the whole-genome transcriptomic profiling, and the remaining sections describe various methods relating to the genomic analyses of the DGRP data.

Genomic feature sets:

We defined genomic feature sets based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, protein–protein interaction networks, and gene expression results from the experiments with the WT males. Using FlyBase annotation “v2016_05” (available at http://flybase.org/) the DGRP SNPs were linked to known genes (within the open reading frame). Genes were then aggregated into KEGG pathways using the BioConductor (Gentleman et al. 2004) package org.Dm.eg.db (Carlson 2015). We obtained known protein–protein interactions from the STRING database (von Mering et al. 2005; Szklarczyk et al. 2015) and restricted the networks to those having a combined score > 0.699 [the score is a combined probability from the different sources of interaction evidence corrected for the probability of observing a random interaction (von Mering et al. 2005)]. Protein identifiers (IDs) were converted to gene symbols using the org.Dm.eg.db package (Carlson 2015) for BioConductor (Gentleman et al. 2004). SNPs were aggregated based on the genes within the gene networks (a total of 7472 subnetworks encompassing 650,766 segregating SNPs). Because we created the subnetworks of directly interacting genes, there is some redundancy, such that genes in one network can also occur in other gene networks. However, this redundancy will be helpful in identifying the particular subset of genes (gene networks) harboring causal SNPs. Finally, we created feature sets based on the WT gene expression profiles by aggregating transcripts according to their P-values (P < 0.001, P < 0.005, P < 0.01, P < 0 0.05, P < 0.1, P < 0.2, and P < 0.5), and whether the transcripts were up- or downregulated, and then linked those transcripts to the SNPs in the DGRP lines.

Gene set enrichment analyses:

KEGG pathways and gene networks were tested for enrichment of transcripts with extreme expression profiles. We used a standard enrichment analysis implemented in the qgg package (Rohde et al. 2019), where for each gene set we computed a summary statistic, , where is the number of transcripts within each gene set and represents the ith t-test statistic from the gene expression analysis (Ackermann and Strimmer 2009; Maciejewski 2014). Then, for each gene set, the summary statistic was compared to a random set of transcripts to assess statistical significance. For each gene set, an empirical P-value was computed as the probability that the observed summary statistic was larger than 10,000 summary statistics computed based on a gene set that was generated by randomly sampling the same number of transcripts as the true gene set. Finally, P-values were adjusted for multiple testing using a false discovery rate < 0.05.

Estimating quantitative genomic parameters:

To determine the existence of genotype-by-treatment interaction in the DGRP, we fitted a bivariate linear mixed model (Schaeffer 1984; Mrode 2005),(3)where the notation is similar to Equation 1, except the subscripts 1 and 2 denote the control and MPH treatment, respectively. The random genetic effects, , and residuals, , were assumed independent and distributed as,(4)Genotype-by-treatment interaction can be assessed in terms of the genetic correlation between the traits (Lynch and Walsh 1998), which can be obtained from the estimated variance components from the bivariate model as (Mrode 2005). If is significantly different from 1 (determined as + SE 1.645 < 1 for P < 0.05), then there is support for a genotype-by-treatment interaction (Lynch and Walsh 1998).

Locomotor response to MPH treatment was quantified as the within-genotype difference in adjusted phenotypic values: , where is the mean of the within-genotype adjusted phenotype for the control and MPH treatment.

The proportion of phenotypic variation of the locomotor response to treatment explained by SNP variation (at MAF > 0.05) was estimated as , where the variance components were estimated by fitting the following,

.(5)

Genomic prediction models:

To predict unobserved phenotypes from observed genotypes, a certain proportion of the data were masked (i.e., became unobserved). A model was then trained on the unmasked data (training data, t ) and the estimated parameters were used to predict the phenotypes in the masked data (validation data, v ). We trained the model on 90% of the data and made predictions using the remaining 10%; this procedure was performed in 50 different random data subdivisions.

First, a null model [genomic best linear unbiased prediction (GBLUP), similar to Equation 5] was fitted containing one random genetic effect for the training data,(6)and the predicted genetic effects in the validation set were computed as,(7)Then, for all feature sets, a GFBLUP (genomic feature best linear unbiased prediction) model was fitted,(8)containing two random genetic components: f and r. The term f is the feature, which captures the genetic variance of the genomic feature, and r is the remaining, capturing the genetic variance not captured by the feature group. The random genetic effects were defined as and , where and . For the GFBLUP model, the total genetic effects in the validation were computed as,(9)The performance of each model was quantified as the correlation between observed and predicted genetic values (here, this corresponds to the adjusted line means); thus, the predictive ability (PA) was . To identify genomic features that significantly increased the PA, each GFBLUP model was compared to the GBLUP model using Welch’s t-test (Welch 1947). Subsequently, all Welch’s t-test P-values were adjusted for multiple testing (false discovery rate) and the significance threshold was set to Padj < 0.05.

The proportion of total genetic variance explained by the genomic feature was estimated as: . This measure is highly dependent on the ability of the model to truly separate the variance between f and r. For small data sets, separating the variance between f and r can be challenging, which can lead to overestimation of . An alternative measure is the proportion of explained variance (PEV), which can be obtained from the relationship between the maximum PA and heritability; the maximum PA is (Mrode 2005; Goddard 2009). Here, we use and the PEV was then estimated as , where was estimated from the parameters obtained in Equation 5, and PA was obtained either from the GBLUP (Equation 6) or the GFBLUP model (Equation 8).

Partitioning of predictive gene sets:

For the gene networks, the GFBLUP model will likely identify a large set of genes that collectively contribute to improved predictive performance. However, it is unlikely that all of those genes contribute equally to the predictive performance. Therefore, we used the covariance association test (CVAT) (Rohde et al. 2016a). This approach allowed us to rank and restrict the combined set of genes to those genes that provide the strongest contribution to the predictive performance (Rohde et al. 2018). The CVAT considers the covariance between the genetic effects of the entire predictive gene network and the genetic effects at the gene level ,(10)where is the genetic effects from a GFBLUP model (Equation 8), where the feature set contains the entire set of genes from the predictive gene networks, and is computed as . The ranking of the genes was based on empirical P-values, by obtaining an empirical distribution of using a circular permutation approach (Cabrera et al. 2012). To retain the correlation structure among the SNPs, the genome was considered circular, and in each iteration of the permutation the genome was rotated such that each SNP would receive a shifted SNP effect. A new, permuted was computed, which this was repeated 10,000 times, and the P-values were then derived from a one-tailed test of the proportion of permuted values that were larger than the observed (Rohde et al. 2016a, 2018; Sørensen et al. 2017).

Data availability

The DGRP genotypes, information on Wolbachia infection status, and chromosomal inversions can be accessed via the website http://dgrp2.gnets.ncsu.edu. The raw and RMA-normalized gene expression data of WT flies are available at the Gene Expression Omnibus under accession number GSE121643. The phenotypic data of the DGRP experiments, the WT experiments, and the gene expression knockdown experiments are available in File S1, File S2, and File S3. Supplemental material available at Figshare: https://doi.org/10.25386/genetics.9368636.

Locomotor and transcriptomic effects of MPH in a WT D. melanogaster population. (A) Diagram of the experimental setup used to quantify the effect of MPH treatment on locomotor activity. Flies were exposed to SUC or MPH, and after 24 hr on the capillary feeding assay, individual flies were moved to the activity chambers where individual movement tracks were obtained using video tracking. (B) Comparison of activities of WT flies subjected to SUC (blue) or MPH treatment (red). Sample sizes (N) and the ANOVA P-value are shown. Scatter points are the observed activities adjusted for experimental factors. (C) Gene expression analysis comparing adult males exposed to SUC or MPH. The volcano plot shows the Benjamini–Hochberg-adjusted P-values (Padj) as a function of fold change in gene expression. Each point represents one transcript, where points highlighted in purple are transcripts that were significantly downregulated and points in orange are transcripts that were significantly upregulated in flies exposed to MPH, relative to flies exposed to SUC. Gene symbols are shown for the 10 transcripts with the highest and lowest fold changes (among the significantly differently expressed genes; full list is available in Table S5). Significance level was set to Padj < 0.001 (horizontal red line) (D) Gene set enrichment analysis of KEGG pathways. Only significant enrichments are shown (Padj < 0.05, Table S6 for all results from the enrichment analysis). Color-coding represents the KEGG pathway classes. KEGG, Kyoto Encyclopedia of Genes and Genomes; MPH, methylphenidate; SUC, sucrose; WT, wild-type.

Genotype-specific response to MPH

We quantified the locomotor activity of 172 DGRP lines exposed to MPH or the control treatment (Table S7), and, in contrast to the WT mass population, we found no overall effect of MPH on activity (P = 0.99, Figure S4A in File S1). However, the genetic correlation of activity between the two treatments was significantly different from one suggesting the presence of a genotype-by-treatment interaction (Lynch and Walsh 1998). On average, flies exposed to SUC ingested a larger volume than flies exposed to MPH (Figure S3 in File S1), but this was uncorrelated with their locomotor activity pattern (Figure S3 and Table S3 in File S1). In addition, the estimated quantitative genetic parameters for DGRP activity at the two treatments did not give rise to concerns related to differential ingestion in the MPH and control flies, since the estimated variance components were very similar, particularly the within-line variances (Table S4 in File S1), suggesting that the proportion of phenotypic variation not accounted for was equal in the two treatment groups.

We quantified the locomotor behavioral response to MPH as the within-DGRP-line difference between the activities of flies exposed to the control and MPH treatments (Table S7). The proportion of phenotypic variation in behavioral response to MPH explained by common SNPs was estimated as (Table S4 in File S1). That is, we found evidence for natural genetic variation for the locomotor response to MPH treatment (Figure 2A). Pearson’s correlation coefficient between the locomotor activity of flies exposed to the control treatment and the locomotor behavioral response to MPH was significant, and positive (P-value < 0.001, Figure S5 in File S1), suggesting that the most active DGRP lines also have the strongest behavioral responses to MPH and thus become less active by MPH treatment, as anticipated from the known effects of MPH in humans.

Genotypic characterization of response to MPH. (A) The locomotor activity response to MPH exposure is dependent on the DGRP genotype. The ranked distribution of locomotor response to treatment for 172 DGRP lines is shown. The purple area shows the DGRP lines that became less active upon MPH treatment, whereas the green area contains the DGRP lines that became more active. (B) Prediction of DGRP locomotor activity after MPH treatment using gene sets defined by levels of differential expression of genes in the WT mass population after exposure to MPH. Color-coding signifies if the sets of genes were up- (orange) or downregulated (purple) in the WT males. Asterisks indicate significantly improved prediction compared to the null model using all available SNPs (gray bar). (C) Effect of gene expression knockdown on locomotor activity. Each bar represents the locomotor response to MPH treatment of the 14 UAS-GAL4 knockdown lines and the two corresponding control lines (gray bars; only the gene DPY-30L1 had the GD control line). Error bars are the SE of the difference of the means. Asterisks denote statistical significance differences from the respective control line: ns, nonsignificant; * P < 0.05, ** P < 0.01, and *** P < 0.001. Gene names that are underlined are the genes with known functions related to histone-modifying processes (Table S10). DGRP, Drosophila Genetic Reference Panel; MPH, methylphenidate; UAS, upstream activating sequence; WT, wild-type.

To investigate if the set of differentially expressed transcripts identified in the WT mass population was important in forming the DGRP behavioral response to MPH, we used an integrative genomic prediction approach: the GFM (Edwards et al. 2016; Sarup et al. 2016; Fang et al. 2017; Rohde et al. 2017, 2018, 2019). This enabled us to predict activity levels after MPH treatment using the genomic data of the DGRP (Figure S1B in File S1). In short, the GFM approach combines genetic data with external biological information to estimate the genetic contribution of SNPs located within predefined sets of SNPs, here referring to sets of differentially expressed genes identified in the WT mass population. Thus, we sought to determine whether the SNPs within sets of differentially expressed transcripts could predict activity following treatment with MPH. Transcripts that were upregulated in the WT mass population after MPH exposure (Padj < 0.005 and < 0.01) significantly improved the predictive performance compared to the null model, which did not utilize prior information (Figure 2B). This suggests that the set of genes that was upregulated due to MPH in the WT population was important for explaining variation in locomotor activity after the DGRP lines ingested MPH. Importantly, we identified a set of genes that was responsive to MPH treatment in the WT mass population and determined that the set of differentially expressed genes was predictive for how the DGRP lines responded to MPH.

To understand the genetic mechanisms of the genotype-specific response to MPH treatment in the DGRP (Figure 2A), we defined gene networks based on directly interacting proteins (7472 networks in total) and used these gene sets to predict the DGRP response to MPH using the GFM approach; thus, searching for gene networks that could predict the DGRP behavioral response to MPH (Figure S1B in File S1). We identified 87 gene networks that, compared to the null model, significantly improved the PA of the behavioral response to MPH (Table S8). The 87 predictive gene networks contained 1737 genes in total, which is too large a number of genes for functional validation. Based on two selection criteria (see below, Figure S1B in File S1), we restricted the list to 36 genes. These genes were further investigated by conditional gene expression knockdown using the binary UAS-GAL4 system (Table S1 in File S1). As a first criterion for selecting genes for functional validation, we ranked the genes within the entire predictive gene network according to their relative contribution to the overall genomic variance (using CVAT). This approach ranks the genes according to explained variance within the set of genes and provides P-values for all the genes within the network, indicating if the SNPs within those genes capture significantly more genomic variance than a random set of SNPs of the same size (Figure S1B in File S1 and Table S9). Among genes that captured significantly more genomic variance than random SNP sets (P-values < 0.05), we selected the top 20 genes for functional validation (that were available as UAS-RNAi). The second criterion was based on the observation that ∼20% of the bait genes (i.e., the central gene of each subnetwork, which that particular network was based on) in the 87 predictive gene networks had common biological functions; namely, they were involved in histone-modifying processes (Figure S1B and Table S10 in File S1). After crossing a ubiquitous GAL4 driver to the UAS-RNAi lines, 14 of the 36 crosses produced viable offspring (Table S1 in File S1). We quantified locomotor activity for the 14 UAS-GAL4 lines, including the respective control lines, after exposure to either the control treatment or MPH, and tested if the behavioral responses of the knockdown lines were significantly different from the behavioral responses of the control lines. Ten of the tested knockdown lines showed a significantly altered response to treatment compared to the control lines; escl, CG18418, and Myo61F resulted in a lower activity pattern compared to the controls, whereas knockdown of Yeti, mRpS28, Cfp1, Mnn1, Sod1, Smurf, and Jarid2 produced flies with increased activity compared to the controls (Figure 2C). These findings support the idea that those genes are important in explaining the observed variability in behavioral responses to MPH treatment.

Cross-generational effects of MPH

The fact that histone-modifying genes were commonly represented among the bait genes in the predictive gene networks could imply epigenetic regulation induced by MPH. Because it is known that histone modifications may lead to cross-generational effects (Vassoler and Sadri-Vakii 2014), we investigated if MPH induced cross-generational behavioral effects on locomotor activity. We exposed virgin males and females from the WT mass population to the control treatment, or to MPH, and crossed males and females in all four combinations of treatments (i.e., ♀SUC ×♂SUC, ♀MPH ×♂MPH, ♀SUC ×♂MPH, and ♀MPH ×♂SUC; Figure S1C in File S1). We then quantified the activities of F1 male offspring kept on standard food. Independent of which parent was exposed to MPH, the F1 male offspring had altered behavioral phenotypes compared to the F1 offspring of parents receiving the control treatment (Figure 3A). The strongest effect was observed when the dams were exposed to MPH and the sires to the control treatment. This suggests that exposure to MPH results in some heritable modifications, such as histone modifications, that have behavioral consequences in F1 flies. The average activities of the F1 control flies were higher than the WT control flies (Figure 1B and Figure 3A), likely because the F1 flies were fed regular food and not fed using the CAFE assay.

Cross-generational effects of MPH exposure. (A) Comparison of locomotor activity of untreated F1 males from parental crosses exposed to SUC or MPH in different combinations. The color of the parental individuals indicates the type of treatment. Blue: SUC; Red: MPH. Sample sizes (N) and the ANOVA P-values are shown. Scatter points are the observed activities adjusted for experimental factors. (B) Gene expression analysis comparing F1 offspring from the cross between females and males from the control treatment (♀SUC×♂SUC) to the offspring from the cross between females exposed to MPH and males exposed to the control treatment (♀MPH×♂SUC). The volcano plot shows the Benjamini–Hochberg-adjusted P-values (Padj) as a function of fold change in gene expression. Each point represents one transcript, where points highlighted in purple are transcripts that were significantly downregulated and points in orange are transcripts that were significantly upregulated in F1 males from ♀MPH×♂SUC, compared to F1 males from ♀SUC×♂SUC. Significance level was set to Padj < 0.05 (horizontal red line). (C and D) Prediction of DGRP locomotor activity after SUC and MPH treatment, respectively, using gene sets defined by levels of differential expression of genes in the WT F1 males. Color-coding signifies if the set of genes was up- (orange) or downregulated (purple)in the WT strain. Asterisks indicate significantly improved prediction compared to the null model using all available SNPs (gray bars). DGRP, Drosophila Genetic Reference Panel; MPH, methylphenidate; SUC, sucrose; WT, wild-type.

We conducted whole-genome transcriptomic profiling of WT F1 males, and compared the transcription levels of offspring from the cross between parents that received the control treatment to offspring from the cross where only the dams were exposed to MPH (here we observed the largest numerical difference, Figure 3A). With a less-stringent significance threshold than for the WT males, we only identified two transcripts (Lcp1 and GC16775) that were significantly differentially expressed (Padj < 0.05, Figure 3B and Table S11).

Using gene sets based on the degree of differential expression in the F1 offspring, we significantly improved the prediction of locomotor activity in the DGRP for both treatments (Figure 3, C and D and Figure S1C in File S1), which indicates that transcripts with the most extreme expression are important for locomotor activity. This was also supported by the enrichment of transcripts with the largest differential expression profiles in the circadian rhythm KEGG pathway (ID 04711, Table S6), a component of which is locomotor activity (Klarsfeld et al. 2003; Rosato and Kyriacou 2006).

When comparing expression levels of transcripts in MPH-treated and control flies in the WT male and the WT F1 male populations (Figure S1D in File S1), we identified 10 transcripts with high probabilities of being differentially expressed in both WT males and WT F1 males (Figure 4A), which represents a larger proportion than expected under the null model of independence (P < 0.0001). This indicates a common response to MPH across generations, which becomes clearer when comparing the results from the gene set enrichment analyses of the gene networks, where 106 gene networks were significantly enriched for transcripts with divergent expression profiles in both generations (Figure 4B and Table S12). Of the 87 gene networks that improved the accuracy of predicting the DGRP response to MPH, 8 showed enrichment of transcripts with extreme expression profiles in the WT males and WT F1 males (Figure 4C). Finally, combining the 87 predictive gene networks into one large network revealed enrichment of transcripts with altered expression in both WT males and WT F1 males (Figure 4D), and the subset of histone-modifying genes was only enriched for transcripts with altered gene expression in the WT F1 individuals (Figure 4D).

Cross-generational transcriptomic analyses. (A) Comparison of the Pr(DE) from the comparison of the WT males and the WT F1 males. Each dot represents one transcript, where the color-coding indicates if the transcript has a probability > 0.5 of being differentially expressed: WT F1 (blue), WT (green), and both WT F1 and WT males (orange). Red dashed lines indicate a probability > 50% of being differentially expressed. The proportion of transcripts with a probability > 0.5 of being differentially expressed is significantly larger than expected under the null hypothesis (P < 0.0001). (B) Number of significant gene networks (Padj < 0.05) enriched for transcripts with divergent expression profiles in the WT males and WT F1 males. Of the 7,472 gene networks, 106 were significantly enriched in both WT males and WT F1 males. Using an empirical permutation approach, we found that the overlap of gene networks was highly significant (P < 0.0001). (C) Gene networks that increased the predictive performance in the DGRP that were also enriched (Padj < 0.05) for differentially expressed transcripts in the WT males (blue) and the WT F1 males (green, marginally nonsignificant, Padj = 0.06). (D) Gene set enrichment analysis for differentially expressed transcripts of the full gene network that increased the predictive performance in the DGRP (green), the subnetwork of histone-modifying genes (orange), or the set of RNAi-validated genes (purple) for WT males (empty) and WT F1 males (hatched). DGRP, Drosophila Genetic Reference Panel; PPI, predictive performance indicator; Pr(DE), probability of a transcript being differentially expressed; RNAi, RNA interference; WT, wild-type.

Discussion

In this study, we performed comprehensive, population-scale molecular genetic characterization of the response to the pharmacological drug MPH, the active compound used to alleviate symptoms associated with ADHD (Kimko et al. 1999). Utilizing a natural outbred population (WT) and a highly inbred reference population (DGRP), we investigated the effect of treatment with MPH at the behavioral and molecular level, within and across generations, by applying an integrative genomic approach allowing us to leverage both database and molecular genetic information.

Similar to findings from other animal models (Choi et al. 2014), we found that MPH treatment led to increased mean locomotor activity (Figure 1B) and major changes in the transcriptome (Figure 1C), for example in carbohydrate metabolic pathways (Figure 1D). Transcriptomic changes can occur within minutes or hours, and in response to most environmental changes; thus, it could be that some of the 64 differentially expressed transcripts are not directly related to the MPH response, but more generally to the ingestion of different food sources. However, transcriptomic responses to drugs have been reported in animal models (Yano and Steiner 2005; Adriani et al. 2006; Beverley et al. 2014) and in humans (Schwarz et al. 2015), testing alcohol (Ghezzi et al. 2013) and methamphetamine (Sun et al. 2011) for example. Using a network-based approach, Sun et al. (2011) found that methamphetamine, another commonly used drug in the treatment of ADHD (López 2006; Goodman 2007), led to alterations in the expression of genes involved in carbohydrate metabolism, which is consistent with results from our KEGG pathway analysis (Figure 1D and Table S6). In total, Sun et al. (2011) identified 68 transcripts that were differentially expressed in response to methamphetamine treatment, of which 22 were among the top differentially expressed transcripts identified in our study (Table S13 in File S1). Other studies have found similar responses to insecticides (Helvig et al. 2004; Daborn et al. 2007; Wan et al. 2014), nicotine (Passador-Gurgel et al. 2007), and caffeine (Bhaskara et al. 2006) (Table S13 in File S1), suggesting a general transcriptomic response to xenobiotic compounds, such as environmental toxins, stimulants, and many other environmental stressors.

Response to treatment depends on the genotype

In humans, several studies investigating the effects of MPH have found evidence for an association between the MPH response and specific alleles (Winsberg and Comings 1999; Kirley et al. 2003; Stein et al. 2005; Joober et al. 2007; Purper-Ouakil et al. 2008; McGough et al. 2009). The focus has mostly been on SLC6A3, a dopamine transporter. Here, we used the DGRP to investigate if the behavioral response to MPH was genetically variable in the D. melanogaster model and, further, if the set of differentially expressed transcripts identified in the WT mass population could predict the behavioral response to MPH in another D. melanogaster population. In contrast to what we found in the WT mass population, we did not observe an overall effect of MPH on locomotor activity across DGRP lines (Figure S4A in File S1). This might be explained by the two populations having different population histories, originating from different environments, and one consisting of a number of homozygous lines and the other being more genetically diverse. Evolutionary forces have likely resulted in allele frequency differences in the two populations; therefore, the causative effects of genetic variants could also be different.

Although we did not find a significant primary effect of MPH on locomotor activity in the DGRP system, we found that, following SUC treatment, the most active DGRP genotypes had the strongest behavioral response to MPH and became less active upon MPH treatment (Figure S5 in File S1). In addition, we provide evidence for a genotype-specific response to MPH (Figure 2A). As mentioned previously, studies with a mouse model for ADHD, the SorcS2 knockdown line, showed that when a WT strain was exposed to amphetamine it became hyperactive, whereas the knockdown line showed decreased activity (Glerup et al. 2014). We observed the same pattern: WT flies showed increased locomotor activity when exposed to MPH (Figure 1B), and gene expression knockdown of escl, CG18418, and Myo61F resulted in decreased activity (Figure 2C). That gene expression knockdown of Yeti, mRpS28, Cfp1, Mnn1, Sod1, Smurf, and Jarid2 resulted in increased locomotor activity does not contradict the validity of our results, but further supports our findings because we have shown that the response to MPH is highly genotype-dependent (Figure 2A and Table S4 in File S1).

The set of genes that became upregulated due to MPH was useful for predicting DGRP activity after MPH exposure (Figure 2B). This was evident when we assessed DGRP activity after MPH treatment, where the set of upregulated genes identified in the WT mass population significantly increased the predictive performance from 0.12 to > 0.2 (Figure 2B and Table S4 in File S1). Using an average gene expression difference to define sets of genes in our feature model might not fully capture the genes that are driving the differences in response among the lines given expected variation in expression levels within lines (Huang et al. 2014). However, by using the average we actually employ a conservative approach, and given that we get a positive result (Figure 2B), we argue that this result and the approach used elaborate our understanding of the effects of MPH across several populations.

We identified 87 highly interconnected networks, which significantly improved the PA of the behavioral response to MPH (Table S8). Each network consists of a bait gene, the gene that all genes within that network interact with, and 20% of the bait genes within the predictive networks are known to be involved in histone-modifying processes (Table S10 in File S1), including epigenetic regulation. Epigenetic mechanisms have previously been linked to neuropsychiatric disorders, including ADHD (Archer et al. 2011), and from studies with animal models it has been shown that different drugs (e.g., cocaine, sedatives, and psychoactive drugs) impose histone modifications (Black et al. 2006; Wang et al. 2007; Bingsohn et al. 2016). Ubiquitous gene expression knockdown of six histone-modifying candidate genes supported the involvement of those genes in the response to MPH (four were significant, Figure 2C). In total, we successfully validated 71% (10 out of 14) of the investigated putative candidate genes (Figure 2C), providing functional evidence for the involvement of those genes in explaining the observed variability in behavioral responses to MPH.

MPH has cross-generational behavioral effects

Cross-generational effects can be paternal, maternal, or a combination of both. They can also be of an epigenetic and transient nature, or they can be epigenetic and heritable, or some combination thereof (Kristensen et al. 2019). We currently lack the data to pinpoint the relative contribution of these mechanisms to the observed cross-generational effects. However, our findings suggested that histone-modifying genes contributed to the variability in the DGRP genotype-specific response to MPH. Thus, to investigate if MPH induced nongenetic changes that could be inherited across generations, we followed the phenotypic response to MPH from WT parental lines to WT F1 offspring. Independent of whether both WT parental sexes or only one parental sex were exposed to MPH, the untreated WT F1 offspring showed increased activity compared to offspring from a cross between WT parents from the control treatment (Figure 3A). This is consistent with findings from other studies using animal models to investigate cross-generational effects of psychoactive drugs (Bingsohn et al. 2016), alcohol (Lee et al. 2013; Yohn et al. 2015), environmental toxins (Carvan et al. 2017; Knecht et al. 2017), and drug stimulants (Yohn et al. 2015).

Although the differences in gene expression of the WT F1 males from the cross between SUC-treated parents and the female parents exposed to MPH, and male parents exposed to SUC (♀MPH×♂SUC), were less pronounced than for the WT males, the transcripts with the numerically lowest P-values (here downregulated) increased the PA for DGRP locomotor activity before (from 0.14 to 0.23) and after (from 0.12 to 0.25) MPH exposure (Figure 3, C and D). Thus, despite the fact that the differences in gene expression of WT F1 males were less pronounced than for the WT males, the genes that do alter expression levels are important for locomotor activity, and the transcriptional response seems to share characteristics between WT males exposed to MPH and untreated WT F1 males originating from parents exposed to MPH. This is based on the observations that: (1) a common set of 10 genes had a high probability of being differentially expressed (> 50%) in WT males and WT F1 males (Figure 4A); (2) a large number of gene networks were enriched for transcripts with extreme expression profiles in both the WT males and the WT F1 males (Figure 4B); and (3) the combined gene network, comprised of the 87 individual networks, increased the PA in the DGRP, and was enriched for transcripts with extreme expression profiles in both WT males and WT F1 males (Figure 4D). The set of histone-modifying genes was enriched for transcripts with altered expression in the WT F1 males and not the WT males (Figure 4D). This latter result further supports the suggested cross-generational phenotypic effects through histone-modifying processes. This is an interesting observation, supporting recent findings that epigenetic mechanisms are linked to ADHD (Archer et al. 2011) and that drugs can impose histone modifications (Black et al. 2006; Wang et al. 2007; Bingsohn et al. 2016). Thus, the importance of and background for the epigenetic inheritance of responses to drug exposure, and the expression of mental disorders in general, should be investigated in much further detail in future studies.

Integrative genomics to uncover the genetic landscape of treatment response

Understanding the genetic basis for complex traits and diseases has a long history, but discoveries based on information obtained from multiple layers of biological organization, utilizing genomic and computational approaches, are still in their infancy. In this study, we provide an improved understanding of the genetic basis underlying variability in response to the psychostimulant MPH utilizing Drosophila and integrative genomic approaches. Leveraging knowledge on gene interactions and transcriptomic profiles, we identified candidate genes that, when gene expression was disrupted, altered the behavioral response to MPH. We further increased the accuracy of predicting the behavioral response to MPH from genotype data by incorporating gene interactions and gene expression profiles. These findings and the approach taken provide results, and a pipeline that generates hypotheses, that can be tested and validated in other species (e.g., humans) to elucidate the genetic mechanisms involved in genotype-specific responses to pharmaceutical intervention. We fully acknowledge the limitations of linking our results directly to humans and also the multiple shortcomings related to, for example, the generality of our results given the strong genotype-specific responses, the impact of differential intake of food from the contrasting treatments (control and MPH), and the fact that only one dose and exposure time were investigated. Still, we propose that the methodological approaches taken are of general interest for studies of the genetic architecture of complex traits, and that the results enlighten our knowledge about responses, including transgenerational ones, to MPH treatment and caution that providing medicine to patients is complex due to highly heterogeneous responses.

Acknowledgments

The DGRP lines and the tubulin-GAL4 driver were obtained from the Bloomington Drosophila Stock Center [National Institutes of Health (NIH) P40 OD-018537, http://flystocks.bio.indiana.edu) and the UAS-RNAi lines were obtained from the Vienna Drosophila Resource Center (www.vdrc.at). We thank Helle Blendstrup and Susan Marie Hansen for technical assistance, and Field Station Mols, Aarhus University for providing facilities where part of the work was conducted. This manuscript was edited for English Language by Charlesworth Author Services (www.cwauthors.com). This study was funded by grants from the Lundbeck Foundation (R155-2014-1724), the Faculty of Health Sciences and Centre for Integrative Sequencing at Aarhus University, the Danish Strategic Research Council (GenSAP: Centre for Genomic Selection in Animals and Plants, contract no. 12-132452), and the Carlsberg Foundation (2013_01_0949) to D.D., and the Danish Natural Science Research Council through a Sapere aude stipend to T.N.K. (DFF-4002-00036).

, 2016The model beetle Tribolium castaneum can be used as an early warning system for transgenerational epigenetic side effects caused by pharmaceuticals. Comp. Biochem. Physiol. C Toxicol.Pharmacol.185–186: 57–64.doi:10.1016/j.cbpc.2016.03.002

The Genetics Society of America (GSA), founded in 1931, is the professional membership organization for scientific researchers and educators in the field of genetics. Our members work to advance knowledge in the basic mechanisms of inheritance, from the molecular to the population level.