Resistance to peroxide food within the DGRP
We found substantial variation among 180 DGRP lines for survival on H2O2 food. We used a mixed model to study the variation in lifespan and, in addition to significant effects of genotype (random effect, likelihood ratio 2 = 36.2, df = 1, P = 1.7x10-9, Methods), we found that weight was a significant predictor, with larger flies surviving longer on H2O2 (fixed effect, β = 0.49, P < 10-6). After removing the effect of weight, we estimated a broad sense heritability for mean survival on H2O2 food of H2 = 51.1 ±8.9% (mean ±SD).
We compared our measures of H2O2 resistance in the DGRP with traits measured in the DGRP in previous studies. These studies have measured either survival or behavioral responses of the DGRP to two other oxidative stressors, paraquat and menadione [29, 30]. We observed Pearson correlations of r = 0.34 (n = 180, P < 10-4) and r = 0.35 (n = 180 P < 10-4) between the mean survival time reported here and those measured on food supplemented with paraquat or menadione, respectively (Fig S1 [29]). We found no correlation between our measures of H2O2 survival and two behavioral traits, the startle response and climbing, measured by Jordan et al. [30] following chronic (13-16 day) exposure to menadione food (data not shown). We also found that the survival times of the DGRP on H2O2 food correlate highly with survival measured under starvation in two different labs (Fig S1 [26, 38]) as well as data from our own lab (Fig 1a). It is notable that the correlation between H2O2 resistance and starvation was greater than the correlation between H2O2 resistance and any other trait published for the DGRP, including survival on food containing paraquat or menadione (Fig S1). These results suggest that H2O2 food may affect feeding or nutrient assimilation in Drosophila.
To test the possibility that H2O2 affects feeding behavior, we measured the amount of food consumed by flies in the CAFE and dye incorporation assays described in the Methods section [39, 40]. Peroxide reduced feeding in each of three different genotypes tested, including strains with both relatively long and short survival time on H2O2 food (Fig 1b and c). During the 24h feeding period in the CAFE assay, flies exposed to liquid H2O2 food consumed no more than the volume lost due to evaporation in chambers without flies, suggesting that the flies consumed very little H2O2 food over 24h (Fig 1c). This finding suggests that mortality in flies exposed to H2O2 is due at least in part to starvation, in addition to any oxidative stress caused by H2O2 exposure.
Genetic associations with peroxide resistance.
We used a linear regression model in PLINK [41] to test for associations of H2O2 resistance with approximately 1.9 million SNPs with a minor allele frequency (MAF) ≥ 5% and < 30% missing genotypes, while accounting for population structure and a significant effect of the major inversion In(2L)t (Methods). We used the q value approach [42] to control the false discovery rate (FDR) and at 20% FDR, 14 variants (all were SNPs) were associated with resistance (Fig 2a, Table S1). Pairwise linkage disequilibrium (LD) among the 14 SNPs suggests that they associate with H2O2 resistance as seven loci, or groups of SNPs in LD (r2 > 0.5, Fig 2b). With only 180 genotypes, we lack the power to analyze SNP-SNP interactions among these loci. However, our data indicate that H2O2 resistance is polygenic within the DGRP.
To investigate these genetic associations further, we performed gene ontology (GO) analysis, looking for biological processes and signaling pathways that are over-represented among the genes associated with survival on H2O2 food. To identify gene-level associations, many studies use the minimum P-value of all variants in a gene (Pmin). One might expect a bias, such that genes with more variants are more likely to have a smaller (more significant) Pmin by chance alone. Indeed, we found that -log10(Pmin) was positively associated with the number of variants per gene (Fig S2), potentially biasing gene-trait associations in favor of genes with more variants [43]. To test for this bias, we compared the top 200 genes ranked by Pmin with the top 200 genes from ten GWAS of randomly permuted phenotypes. Out of 15,322 gene models, the null expectation for such intersection would be 2.6 genes. In contrast, we found an average of 11.6 ± 2.0 genes (mean ± se) in common across the permutations (2 = 28.3, df = 3, P < 1.1x10-7), consistent with a bias caused by SNP density. To correct for this bias, we used a permutation approach to derive gene-level P-values, Pgene, while also accounting for population structure (Methods). Unlike Pmin, Pgene did not associate with the number of variants per gene, and it reduced the number of false-positives when compared to top genes from GWAS of randomized phenotypes (Fig S2, 2 = 2.2x10-26, df = 3, P = 1). Thus, Pgene increases the accuracy of gene-trait associations. Among the 29 genes at Pgene FDR ≤ 0.5 (Table 1) we used gene ontology (GO) enrichment analysis to identify several enriched biological processes and two enriched biological pathways (FDR ≤ 0.05, Table 2). Most of the enriched processes are nested in hierarchical categories and thus are not independent. Also, the enrichment of 7 of the 9 biological processes, and the endothelin signaling pathway, is due entirely to three genes, each encoding adenylyl cyclase (ACXA, ACXB and ACXE, Table 2). Adenylyl cyclase is involved in several signaling pathways, including G-protein coupled and calcium-based signaling. Separately, the platelet-derived growth factor (PDGF) signaling pathway is enriched (FDR = 0.02) due to three other genes in our dataset, Rab2, Ets21C and the c-Myc-binding protein homolog CG17202.
[Please place Table 1 here]
Table 2. Process and Pathway Enrichment Among Candidate Genes
PANTHER GO-Slim Biological Process
|
Gene Content
|
Input
|
Expected
|
Fold Enrichment
|
P-value
|
FDR
|
activation of adenylate cyclase activity (GO:0007190)
|
32
|
3*
|
0.07
|
46.09
|
4.69x10 -5
|
0.0353
|
adenylate cyclase-activating G-protein coupled receptor signaling pathway (GO:0007189)
|
32
|
3*
|
0.07
|
46.09
|
4.69x10 -5
|
0.0235
|
cAMP-mediated signaling (GO:0019933)
|
41
|
3*
|
0.08
|
35.98
|
9.38x10 -5
|
0.0282
|
regulation of cAMP-mediated signaling (GO:0043949)
|
41
|
3*
|
0.08
|
35.98
|
9.38x10 -5
|
0.0235
|
regulation of adenylate cyclase activity (GO:0045761)
|
42
|
3*
|
0.09
|
35.12
|
1.00x10 -4
|
0.0216
|
regulation of lyase activity (GO:0051339)
|
43
|
3*
|
0.09
|
34.3
|
1.07x10 -4
|
0.0201
|
cyclic-nucleotide-mediated signaling (GO:0019935)
|
50
|
3*
|
0.1
|
29.5
|
1.64x10 -4
|
0.0274
|
regulation of biological process (GO:0050789)
|
1359
|
11
|
2.76
|
3.98
|
3.87x10 -5
|
0.0582
|
biological regulation (GO:0065007)
|
1492
|
11
|
3.03
|
3.62
|
9.15x10 -5
|
0.0343
|
|
|
|
|
|
|
|
PANTHER Pathways
|
Gene Content
|
Input
|
Expected
|
Fold Enrichment
|
P-value
|
FDR
|
PDGF signaling pathway (P00047)
|
47
|
3
|
0.1
|
31.38
|
1.38x10-4
|
0.0213
|
Endothelin signaling pathway (P00019)
|
74
|
3*
|
0.15
|
19.93
|
0.0495
|
0.0256
|
Table 2. Enrichment analysis for the 29 genes associated with peroxide resistance. Gene Content is the number of genes in the respective process/pathway the Drosophila genome (n = 13,767 total gene models). Asterisks indicate categories enriched due to the adenylyl cyclase X gene cluster.
We validated these gene-trait associations by using RNAi to manipulate the expression of six of the 29 candidate genes, nAChRbeta3, Ets21C, ush, Nha1, Jon25Bi and Marcal1 (Pgene < 0.0003 in all cases), and testing the effect of mutation in a seventh candidate NPF (Pgene < 0.0006). Several of the candidates reside near a cluster of six trait-associated SNPs within a 17.5 kb interval on chromosome 2 (Fig 2a). This interval spans several genes, including u-shaped (ush, Pgene = 7.1x10-5), which contains an intronic C/T SNP associated with H2O2 resistance (P = 2.49x10-8) and several other SNPs in LD with this SNP (Fig 2b). None of these SNPs are predicted to alter amino acid sequence or splicing, but instead may affect ush expression. ush has roles in development and growth, including as a negative regulator of PI3K activity within the IlS pathway in the fat body [44]. To validate the effects of ush on H2O2 resistance, we used the RU486-inducible GAL4 GeneSwitch driver S106 to drive RNAi targeting ush in the adult fat body [45]. RU486 treatment of flies carrying both S106 and UAS-RNAi targeting ush resulted in shorter survival times on H2O2 food than the same genotype without RU486 (Fig 3). We saw the same result with two independent RNAi lines, each targeting a different portion of the ush mRNA (Fig 3). We saw no effect of RU486 on the survival of F1 flies carrying the driver and the empty control P-element in the same genetic background as the UAS-RNAi flies, nor an effect of RNAi on lifespan on food lacking H2O2 (see Methods). Knocking down ush in the nervous system with the elav-GAL4 gene switch driver did not affect H2O2 resistance (data not shown), and ubiquitous knockdown of ush appears to be lethal as we failed to recover Act-GAL4/ush-RNAi flies in crosses of the ush-RNAi construct to the constitutive Act-GAL4 driver. We also detected a strong interaction between RU486 treatment and an Act-GAL GeneSwitch driver in our H2O2 food assay and so were unable to test the effect of ubiquitous knockdown of ush in adult flies (data not shown). Similar RNAi of five other candidates did not affect H2O2 resistance (data not shown).
Another candidate, NPF (Pgene = 5.8x10-4) encodes the ligand neuropeptide F (NPF), which controls feeding, ethanol sensitivity and other behaviors [46]. The npfSK1 deletion allele reduced H2O2 resistance when compared to wildtype control flies and did not affect survival on control food (Fig 3 and data not shown). These data demonstrate that manipulation of candidate genes from our GWAS affects the H2O2 resistance phenotype.
Metabolite profiles associated with peroxide resistance.
To investigate the effect of H2O2 on the fly metabolome, and the potential for the metabolome to explain genetic variation in resistance to H2O2, we measured untargeted LC-MS profiles in three biological replicates for each of eight resistant (mean survival time = 106.9h, range = 90.4 - 119.7h) and eight sensitive (mean survival time = 58.9h, range = 53.2 - 70.0h) lines, chosen such that the two groups did not differ substantially in size (Fig S3). These lines were subjected to another survival assay and, 24h after being exposed to H2O2 or control food, samples of flies from each line and treatment were flash frozen for aqueous metabolite extraction, while survival measurements were conducted on the remaining flies.
Global metabolite profiling identified a total of 3028 and 2921 features from positive and negative ionization modes, respectively. Profiles from positive and negative modes were analyzed separately. Principal component analysis (PCA) was applied to metabolomic data (see Methods). For the negative mode data, principal component one (PC1neg) explained ~12.9% of the variation and separated resistant from sensitive flies but did not separate flies based on treatment alone (Fig 4). PCA using positive mode data gave similar results (data not shown). Analysis along PC1 also separates H2O2-treated from control flies among the sensitive genotypes, but not among the resistant genotypes (Fig 4 and data not shown). Principal component analysis thus detected between-group variation in metabolite profiles from sensitive and resistant flies, and further distinguished the effect of treatment on sensitive but not resistant flies.
The separation of resistant and sensitive lines by metabolome profile is striking. However, the genotypes chosen for metabolite profiling were among the extremes of resistance to H2O2. This design raises the possibility that resistant and sensitive lines have a genetic signature of resistance. To determine if genotype could also separate lines chosen deliberately with extreme phenotypes, we carried out PCA and hierarchical clustering on the same lines using their genotypes. We analyzed the first ten genotype PCs, which together account for 68.7% of the variance of >2 million genetic variants in the 16 lines used for metabolomics. These PCs failed to clearly separate the resistant and sensitive flies (Fig 5). Similarly, clustering of 6846 LD-pruned variants from these 16 lines also failed to separate the two phenotypes. For comparison, clustering of 2417 negative mode metabolite features separated 14 of the 16 genotypes into two clusters composed of resistant and sensitive lines (Fig 5).
Metabolic pathways associated with peroxide resistance
We next sought to identify individual features and metabolic pathways that predicted the binary trait of resistance or sensitivity to H2O2 food. Accordingly, we used logistic regression, identifying a large number of metabolite features from flies on H2O2 food or control food whose abundance was associated with trait. Among the features detected in negative mode LC-MS of untreated (control) flies, 252 features associated with trait (FDR < 0.1), 228 of which were also found in the larger set of 637 associated features identified in the metabolome of H2O2 treated flies (Fig S4a). In positive mode data, 115 features associated with resistance in H2O2 treated flies, but no features did so from flies on control food.
Principal component analysis suggested that the metabolome of sensitive lines is altered by H2O2 food, while in resistant lines, the metabolome does not shift in response to H2O2 treatment (Fig 4). To identify the specific metabolites whose abundance was affected by treatment in a trait-dependent way we ran a linear regression model, predicting metabolite level in response to trait, treatment and the interaction between the two. Fig S4b shows the clustering of 138 features that show significant interaction term at FDR < 0.1. These data come from analysis of profiles from negative mode only, as no features from positive mode reached our FDR cutoff. Clustering the z-scores for each feature across samples revealed several groups with similar patterns, including one large group of features in which the sensitive flies, but not the resistant flies, showed a decrease in feature abundance on H2O2 compared to control treatment (Fig S4). This pattern is similar to the separation of samples across the latent variables revealed in PCA, where the metabolome of sensitive flies, but not resistant flies, was affected by treatment (Fig 4).
We used the mummichog software package [47] to identify metabolomic pathways enriched among the features associated with H2O2 resistance, or among features with significant treatment by trait interaction effects. This analysis identified metabolites that were enriched in several pathways, including carbohydrate metabolism, amino acids and their biosynthesis, and folate metabolism (see Additional File 1). Many of these pathways share overlapping metabolites and some pathways are nested within other pathways. We chose a subset of the identified pathways for further analysis based on a variety of criteria. One criterion was the significance of each pathway across the different measures of association. Glycogen and folate metabolic pathways were enriched among the features that were significant in at least three of the four analyses (see Additional File 1). Another criterion was the strength of the identification; we gave higher priority to those features that were uniquely assigned to a particular metabolite or pathway, rather than being ambiguously associated with more than one metabolite, or with a metabolite that is present in several pathways.
Glycogen Metabolism
Among the features that showed significant trait by treatment interaction were those with glycogen metabolism, particularly glucose and the maltodextrins containing from two to six glucose moieties (Fig S5). The next maltodextrin in this series, maltoheptaose, has a mass (1153 Da) outside the range of features measured in this experiment (62 to 1000 Da, see Methods). Most of the maltodextrins show the same general relationship with treatment and trait; their abundance is lower in sensitive flies on H2O2 food compared to unexposed sensitive flies, and their abundance is not affected by treatment in resistant flies (Fig S4b). In contrast to the unique assignment of features to the maltodextrins, the other metabolites of the glycogen pathway, glucose and the disaccharide maltose, had features that could be assigned to other metabolites of identical mass. Together these data suggest that glycogen metabolism differentiates sensitive verses resistant flies.
Resistant flies appear to maintain a glycogen pool after transfer to H2O2 food whereas glycogen is relatively depleted in sensitive flies, suggesting that their survival may be limited by their glycogen reserve (Fig 6b). Metabolite levels were measured 24h after exposure to H2O2, and while the glycogen pool of resistant flies appears to persist until this point, we do not know the kinetics of glycogen levels over the course of their typical survival on H2O2, which averages 107h. Consistent with a role for glycogen in resistance to H2O2 food, flies fed supplemental maltose prior to being exposed to H2O2 showed increased resistance to H2O2 food in a dose-dependent manner (Fig 6). Maltose could increase survival by providing glucose, or perhaps by some other effect as a disaccharide. Supplementing fly diet with glucose but not lactose extends survival time similar to supplemental maltose, which is consistent with the former hypothesis (Fig 6).
The fat body is a site of glycogen storage in Drosophila, and RNAi targeting glycogenin (CG44244), the gene encoding the protein core of glycogen, in the fat body increased the sensitivity of flies to H2O2 food (Fig 7). We found that knocking down glycogenin using the S32 fat-body driver reduced survivorship, while the S106 fat-body driver did not (Fig 7). Glycogen is present in the Drosophila brain, and knocking down glycogenin with the neuron-specific elavGS driver reduced survival on H2O2 substantially (Fig 7) [48].
Folate Metabolism
Metabolites of the folate pathway were overrepresented in features associated with H2O2 resistance (see Additional File 1). The folate pathway is central to the synthesis of several amino acids, nucleotides, secondary metabolites and substrates for secondary modifications (e.g. methylation). We detected features corresponding to metabolites both in the folate pathway as well as in peripheral pathways (e.g. S-adenosylmethionine), suggesting that the activity of the folate pathway differs between sensitive and resistant flies. To investigate the potential role of the folate pathway in H2O2 resistance, we tested the effect of metabolic gene knockdown on survival. Knocking down either CG8665, which encodes 10-formyltetrahydrofolate dehydrogenase (FDH), or pugilist (CG4067), which encodes tetrahydrofolate dehydrogenase (THFDH), in the fat body or in neurons reduced the survival of flies on H2O2 food (Fig 7). In parallel experiments, RU486 pretreatment failed to affect survival in flies carrying the fat body 32GS driver. These data suggest that folate metabolism in the abdominal fat body and neurons is important for survival on H2O2 food.
Supplemental folic acid was shown to reduce the levels of oxidative damage associated with knockdown of parkin in Drosophila neurons [49]. We also tested whether supplemental folic acid would affect survival on H2O2 food and failed to see a significant effect (data not shown).