The Mechanisms Underlying α-Amanitin Resistance in Drosophila melanogaster: A Microarray Analysis

Abstract

The rapid evolution of toxin resistance in animals has important consequences for the ecology of species and our economy. Pesticide resistance in insects has been a subject of intensive study; however, very little is known about how Drosophila species became resistant to natural toxins with ecological relevance, such as α-amanitin that is produced in deadly poisonous mushrooms. Here we performed a microarray study to elucidate the genes, chromosomal loci, molecular functions, biological processes, and cellular components that contribute to the α-amanitin resistance phenotype in Drosophila melanogaster. We suggest that toxin entry blockage through the cuticle, phase I and II detoxification, sequestration in lipid particles, and proteolytic cleavage of α-amanitin contribute in concert to this quantitative trait. We speculate that the resistance to mushroom toxins in D. melanogaster and perhaps in mycophagous Drosophila species has evolved as cross-resistance to pesticides, other xenobiotic substances, or environmental stress factors.

About the Authors:

Chelsea L. Mitchell

Affiliation: Department of Biological Sciences, Michigan Technological University, Houghton, Michigan, United States of America

Michael C. Saul

Affiliation: Department of Zoology, University of Wisconsin-Madison, Madison, Wisconsin, United States of America

Liang Lei

Affiliation: School of Forest Resources and Environmental Science, Michigan Technological University, Houghton, Michigan, United States of America

Hairong Wei

Affiliation: School of Forest Resources and Environmental Science, Michigan Technological University, Houghton, Michigan, United States of America

Thomas Werner

* E-mail: [email protected]

Affiliation: Department of Biological Sciences, Michigan Technological University, Houghton, Michigan, United States of America

Introduction

How species respond to changes in their environment is a central question in biology. Insects and mammals deploy similar genes and detoxification mechanisms to defend against poisons that are present in their prey or in the environment. These include the avoidance of toxic parts of their diet, the excretion, sequestration, metabolic breakdown of the toxins, and mutations in the target proteins to avoid toxin binding [1]. Some of the most striking natural examples of toxin resistance are snake species that feed on poisonous amphibians [2], caterpillars that sequester plant alkaloids in their bodies to deter predators [1], and toxin-resistant soft-shell clams that store algal toxins in their bodies, causing paralytic shellfish poisoning in people who eat the clams [3]. Apart from these natural examples, the use of pesticides against insects has caused very rapidly evolving toxin resistance responses in many pest species [4]–[8], costing the US billions of dollars per year in crop damage and pesticide production [9].

Out of the vast number of eukaryotic organisms that live on our planet, a few dozen of mycophagous Drosophila species are able to breed in a variety of very toxic mushrooms, including the deadly poisonous species Amanita phalloides (Death Cap) and Amanita virosa (Destroying Angel). Among other toxins, these mushrooms contain α-amanitin as their principal toxin, which inhibits the function of RNA-polymerase II and thus brings all mRNA transcription to a halt [10]. These resistant Drosophila species can develop on α-amanitin-containing laboratory food [11], [12], showing that the resistance mechanism is not due to the avoidance of toxic parts of the mushrooms. Furthermore, the RNA-polymerase II of all tested mushroom-feeding Drosophila species is as sensitive to α-amanitin as it is in sensitive Drosophila species [13], showing that target mutations in the RNA polymerase II complex are not likely to confer resistance to mushroom toxins in mycophagous Drosophila species.

The model organism D. melanogaster is a non-mycophagous species; i.e., it does not use mushrooms as a natural diet. Thus, D. melanogaster should not encounter toxins in nature that are solely produced by mushrooms, such as α-amanitin. However, three Asian D. melanogaster strains that were collected in the 1960s in Taiwan (Ama-KTT), India (Ama-MI), and Malaysia (Ama-KLM) were shown to be one order of magnitude more resistant to α-amanitin than the sensitive wild-type strain Oregon-R [14]. In these three Asian strains, the resistance to α-amanitin was mapped to two dominantly acting loci: one situated on the left and one on the right arm of chromosome 3. Eighteen years later, a very similar phenomenon was described in a D. melanogaster stock collected in California. This stock showed an increased resistance to α-amanitin and surprisingly, the resistance was mapped to the seemingly same two loci on chromosome 3, as in the three Asian stocks. Even in the Californian stock, both loci acted in a dominant fashion [15]. The Californian study concluded with the identification of two candidate genes that might confer the resistance phenotype: Multidrug resistance 65 (Mdr65) on the left arm and Protein kinase C98E (Pkc98E) on the right arm of chromosome 3. Because PKC98E can phosphorylate MDR proteins [16] and MDR proteins could potentially lead to the excretion of α-amanitin from cells, the question of how D. melanogaster evolved α-amanitin resistance appeared to be answered. Although the proposed scenario is simple and elegant, no conclusive evidence has been brought forward yet that demonstrates that any gene is required or necessary to confer resistance to α-amanitin. Thus, the genes that confer mushroom toxin resistance in D. melanogaster (and all mycophagous Drosophila species) remain elusive.

In this study, we conducted a whole-genome microarray analysis, using an isochromosome stock for chromosomes 2 and 3 of the original α-amanitin-resistant D. melanogaster stock Ama-KTT from Taiwan. We hypothesized that genes involved in the excretion, metabolic inactivation, and/or sequestration of α-amanitin will be identified in our microarray, which can pinpoint to the mechanisms responsible for the α-amanitin resistance phenotype. To our surprise, neither Mdr genes nor Pkc98E were among the up-regulated candidate genes. Instead, we identified genes of the phase I detoxification gene family Cyp (Cytochrome P450), and the phase II Gst (Glutathione-S-transferase) and Ugt (UDP glucuronosyl transferase) gene families, some of which (Cyp6a2, Cyp12d1-d, and Cyp12d1-p) were several hundred-fold constitutively up-regulated in the α-amanitin-resistant fly stock. In addition, we found evidence for the possible involvement of peptidases, lipid particles, cuticular proteins, the Mayor Royal Jelly Protein homolog Yellow, and Salivary Gland Secretion (Sgs) proteins, which could provide additional protection by cleaving or immobilizing α-amanitin, or by blocking its access to cells. Because D. melanogaster does not feed on mushrooms in nature and α-amanitin is solely found in mushrooms, we speculate that the resistance to α-amanitin has evolved as cross-resistance to pesticides or other environmental factors that the flies encountered before they were collected in Asia 45 years ago.

Results

Experimental Design

In two independent studies, a total of four D. melanogaster stocks from Asia and North America were shown to be resistant to the mushroom toxin α-amanitin [14], [15]. For each of these stocks, QTL mapping data suggested that the resistance was conferred by two dominantly acting loci on chromosome 3. Begun and Whitley identified the genes Mdr65 and Pkc98E as possible candidates, with the notion that the resistance phenotype could be caused by a cis-regulatory change in the Mdr65 gene [15]. In order to identify gene-regulatory changes on a whole-transcriptome scale in α-amanitin-resistant D. melanogaster larvae, we performed a microarray study. As starting material, we used the most resistant of the four described α-amanitin-resistant stocks, Ama-KTT [14]. Because the stock could have become heterozygous for the resistance-conferring loci during the past 45 years after being collected in the wild, we created the isochromosome stock Ama-KTT/M/2, which is isogenic for the Ama-KTT chromosomes 2 and 3. Our dose-response data show that the isochromosome stock Ama-KTT/M/2 (LC50 = 2.16 μg/g of food) is at least as resistant to α-amanitin as the original Ama-KTT stock (LC50 = 1.84 μg/g of food) (Figure 1), indicating that at least the majority of the α-amanitin resistance-conferring genes is located on the major autosomes. The multi-balancer stock that we used for the crosses to create the Ama-KTT/M/2 stock was very sensitive to α-amanitin (LC50 = 0.042 μg/g of food, data not shown).

[Figure omitted. See PDF.]

Figure 1. Ama-KTT/M/2 is not less resistant to α-amanitin than Ama-KTT.

Ten first-instar larvae were placed on each α-amanitin concentration. The dose response curve shows the percentage of hatching flies. Error bars indicate the s.e.m. of three replicates.

https://doi.org/10.1371/journal.pone.0093489.g001

We performed a whole-transcriptome gene expression microarray analysis to test what genes are differentially expressed in 1) a constitutive manner and 2) in response to α-amanitin. The complete set of differentially expressed genes (DEGs) can be found in Table S1. The isochromosome stock Ama-KTT/M/2 (LC50 = 2.16 μg/g of food) was used as the experimental stock and has a 77.1 times higher LC50 to α-amanitin than our sensitive control stock Canton-S (LC50 = 0.028 μg/g of food, data not shown). We compared three groups with each other: 1) Canton-S larvae on non-toxic food, 2) Ama-KTT/M/2 larvae on non-toxic food, and 3) Ama-KTT/M/2 larvae that were continuously raised from the first to the third instar on α-amanitin-containing food (at 1.5 μg/g of food, a concentration that is slightly below the LC50 of Ama-KTT/M/2). Groups 1 and 2 were prepared in 5, and group 3 in six biological replicates, each replicate consisting of ten larvae (Figure 2). We compared the gene expression profiles of fully-grown third-instar larvae that have not started wandering yet. For the data analysis, we focused on well-annotated genes that showed expression changes of at least 2-fold with a corrected p-value of less than 0.05. With the exception of the genome enrichment analysis and the gene CG10226, which is a putative Mdr gene, we generally excluded genes from our analysis that solely have a CG or CR gene annotation number.

[Figure omitted. See PDF.]

Figure 2. Graphical representation of the groups of larvae used for the microarray and qPCR analysis.

Groups 1 and 2 (Canton-S and Ama-KTT/M/2) were not treated with α-amanitin, as symbolized by the yellow color. The larvae of group 3 (Ama-KTT/M/2) were treated with α-amanitin throughout their development, as indicated in red. Groups 1 and 2 were collected in five, and group 3 in six biological replicates (ten larvae in each replicate), as illustrated by the number of tubes and microarray chips.

https://doi.org/10.1371/journal.pone.0093489.g002

Genes Encoding Cytochrome P450s, GSTs, and UGTs Are Differentially Expressed in Ama-KTT/M/2

Assuming that gene-regulatory changes underlie α-amanitin resistance in the Ama-KTT/M/2 isochromosome stock, we expected to identify constitutive gene-expression changes in Ama-KTT/M/2 on non-toxic food, as compared to the sensitive control stock Canton-S on non-toxic food (group 2 versus group 1). We used the Plier normalization/summarization and the DEG methods to analyze our single gene microarray data. As a result, we identified 234 genes that were at least 2-fold significantly constitutively up-regulated in Ama-KTT/M/2 (Table S2). Out of these 234 genes, 20 (8.5%) are Cyp, Gst, and Ugt genes, which are all situated on chromosomes 2 and 3 (Table 1). The three most highly up-regulated genes of this group were Cyp6a2, Cyp12d1-d, and Cyp12d1-p, which were between more than 300- to 197.3-fold constitutively up-regulated in the resistant stock. These three genes are expressed in the larval midgut and Malpighian tubules, which are potential detoxification organs [17]. Interestingly, Cyp6a2 expression profiles are correlated with insecticide resistance [18]–[24], while CYP6A2 metabolizes insecticides in enzyme assays [19], [25]. Cyp12d1 is also associated with insecticide resistance [22], [23], [26]–[31] and stress response [29]–[31]. Overexpression of Cyp12d1 increases insecticide resistance [32], and CYP12D1 from the house fly metabolizes insecticides [33]. The remaining 17 detoxification genes identified in our microarray study were 38.9 - 2.1-fold up-regulated and are presented next in the order from the highest to lowest constitutive up-regulation in Ama-KTT/M/2: Ugt36Bb, Ugt86Dd, GstD5, GstE1, GstE6, GstE5, Ugt36Bc, Cyp6a20, Ugt37c1, Ugt36Ba, Cyp4c3, Ugt37b1,Cyp6w1, Cyp305a1, Cyp49a1, GstD8, and GstE9. Some of these genes are associated various phenotypes: Ugt86Dd and Cyp6w1 (inducibility by the xenobiotic phenobarbital) [23], GstD5 and GstE1 (stress responses) [34], [35], GstE5 (insecticide resistance) [36], Cyp6a20 (aggressive behavior) [37]–[39], and Cyp305a1 (ecdysteroid synthesis and lipid storage regulation) [40].

[Figure omitted. See PDF.]

Table 1. Single gene analysis for Ama-KTT/M/2 versus Canton-S on no toxin (group 2 versus 1).

https://doi.org/10.1371/journal.pone.0093489.t001

We were curious to see if the constitutive up-regulation of detoxification genes is a specific characteristic for the α-amanitin-resistant stock Ama-KTT/M/2 or if there are other detoxification genes that show higher expression levels in Canton-S, as compared to Ama-KTT/M/2. Surprisingly, 15 Cyp and Ugt genes were between 2.1 and 186.8-fold lower expressed in the resistant stock Ama-KTT/M/2 than in Canton-S. From the lowest to highest expression difference, these genes are: Cyp12a4, Cyp304a1, Cyp313a2, Cyp12e1, Cyp6t1, Cyp4ac2, Cyp4s3, Ugt86Dj, Cyp4d2, Cyp6a23, Cyp4ac3, Cyp4p2, Cyp28d1, Cyp4d8, and Cyp6a17 (Table 1). Correlative or functional data exists for Cyp12a4 and Cyp4p2 (insecticide resistance) [21], [41], Cyp304a1 and Cyp4d2 (methanol resistance) [42], and Cyp6a17 (thermosensory behavior) [43].

Genes Encoding Cytochrome P450s and GSTs Are Inducible by α-Amanitin

Our next question was what genes are inducible by α-amanitin in the resistant Ama-KTT/M/2 stock as compared to Ama-KTT/M/2 on non-toxic food (group 3 versus group 2). We found that 143 genes were significantly inducible by α-amanitin (Table S3), eleven of which (7.7%) belong to the Cyp and Gst gene families (Table 2). Cyp316a1 was the strongest inducible Cyp gene (11.8-fold) in the resistant stock Ama-KTT/M/2. However, when we compared resistant Ama-KTT/M/2 on toxic food to sensitive Canton S without toxin (group 3 versus group 1), Cyp316a1 was only 1.9-fold (p = 0.0941, Table S1) more expressed in Ama-KTT/M/2 on toxic food, making the 11.8-fold induction within the Ama-KTT/M/2 stock less convincing. The remaining ten Cyp and Gst that were up-regulated by α-amanitin in the resistant stock were induced between 7.2- and 2.0-fold and are listed in the order from highest to lowest induction: Cyp6d2, Cyp4d8, Cyp28d1, Cyp6t1, GstD3, GstD6, Cyp4d2, GstD9, GstD10, and Cyp4d14. Four of these genes, Cyp4d8, Cyp28d1, Cyp4d2, and Cyp4d14, are expressed in the larval midgut and/or Malpighian tubules, suggesting that they could play a role in the detoxification of xenobiotic compounds [17]. Some genes are associated with various phenotypes: Cyp6d2 (camptothecin resistance) [44], GstD6 (oxidative stress response) [45], and Cyp4d2 (methanol resistance) [42]. Notably, both Cyp316a1 and Cyp4d8 are situated at cytological position 66A2, which is relatively close to region 65A10 to which α-amanitin resistance was QTL-mapped in four independent D. melanogaster stocks in the past [14], [15]. We next asked what Cyp, Gst, and Ugt genes were down-regulated in response to α-amanitin in the resistant stock. As a result, nine genes were 2.1- to 3.8-fold down-regulated in response to α-amanitin, which are presented in the order from lowest to highest down-regulation: Ugt37b1, Cyp4c3, Cyp28d2, Ugt86Dd, Cyp6a23, Cyp9b2, Ugt37c1, and Cyp28a5 (Table 2). Out of these, Ugt86Dd is inducible by the xenobiotic phenobarbital [23] and Cyp28a5 by methanol [42]. Some of the most strongly α-amanitin-inducible genes (>300-fold) were the salivary gland secretion genes Sgs1, Sgs3, Sgs5, Sgs7, and Sgs8 (Table S3). We will speculate about their role later.

[Figure omitted. See PDF.]

Table 2. Single gene analysis for Ama-KTT/M/2 on α-amanitin versus Ama-KTT/M/2 on no toxin (group 3 versus 2).

https://doi.org/10.1371/journal.pone.0093489.t002

Mdr Genes Are Neither Constitutively Up-Regulated nor Inducible in Ama-KTT/M/2

In 1982, QTL mapping data suggested that two loci on chromosome 3 of the Asian Ama-KTT, Ama-MI and Ama-KLM stocks confer resistance to α-amanitin in a dominant fashion [14]. Eighteen years later, a Californian D. melanogaster stock showed α-amanitin resistance that was QTL-mapped to virtually the same two loci on chromosome 3 [15]. It was concluded that Mdr65 and Pkc98E were possible candidate genes for causing the resistance. Furthermore, sequence comparisons between the most and the least resistant Californian stocks pointed out differences in the non-coding regions, but not in the coding regions of Mdr65. Thus, if Mdr65 would confer resistance, the prediction was that a cis-regulatory change in the Mdr65 gene is responsible for the resistance α-amanitin. We thus asked the question if Pkc98E, Mdr65 or any other Mdr genes (CG10226, Mdr49, and Mdr50) were either constitutively up-regulated or inducible by α-amanitin in the Ama-KTT/M/2 stock. Comparing group 2 with group 1, Mdr65 showed a statistically significant but very low (1.2-fold) constitutive up-regulation in Ama-KTT/M/2 (Table 1), while Mdr65 was 1.2-fold down-regulated in response to α-amanitin when group 3 was compared to group 2 (Table 2). CG10226, a predicted Mdr gene that directly flanks the Mdr65 gene on the left arm of chromosome 3, showed a statistically significant 1.4-fold constitutive up-regulation in the resistant Ama-KTT/M/2 stock as compared to Canton-S (Table 1), while this gene was 1.2-fold down-regulated in response to α-amanitin (Table 2). The remaining two Mdr genes of D. melanogaster, Mdr49 and Mdr50, are both situated on the right arm of chromosome 2. Mdr49 showed a mere 1.1-fold constitutive up-regulation in Ama-KTT/M/2 (Table 1), and it is 1.1-fold inducible by α-amanitin (the latter value is statistically insignificant) (Table 2). The observed 1.6-fold constitutive induction of the Mdr50 gene was statistically significant (Table 1), and the same gene was significantly 3.3-fold down-regulated in response to α-amanitin (Table 2). Furthermore, Pkc98 is 1.1 times lower expressed in Ama-KTT/M/2 as compared to Canton-S on no toxin (Table 1), while this gene is 1.2-fold inducible by α-amanitin (both values statistically insignificant) (Table 2). In summary, our data show that Mdr genes and Pkc98E were far less than 2-fold (if at all) up-regulated, neither constitutively nor in response to α-amanitin. Mdr genes are thus not likely to confer the α-amanitin resistance, at least not by increasing Mdr gene expression.

We also specifically analyzed the regulation of two transcription factor genes that are known to play a role in regulating responses to xenobiotic factors. Hr96 encodes a nuclear receptor that is involved in xenobiotic responses in D. melanogaster [46]. Our data in Table 1 show that Hr96 is 1.6-fold constitutively higher expressed in Canton-S (group 1) than in Ama-KTT/M/2 on no toxin (group 2). In response to α-amanitin, the Hr96 gene was 1.2-fold (statistically insignificant) down-regulated in Ama-KTT/M/2 (group 3 versus group 2, Table 2). The other gene of interest was the leucine zipper transcription factor cnc, which is known to activate oxidative stress and detoxification responses in D. melanogaster [47], [48]. In our microarray, the cnc gene is 1.2-fold constitutively higher expressed in Canton-S (group 1) than Ama-KTT/M/2 (group 2, Table 1), while α-amanitin treatment caused a1.2-fold induction in the Ama-KTT/M/2 stock (group 3 versus group 2, Table 2), thus bringing cnc gene expression to the same level that was observed in Canton-S without toxin. These results, at least at the transcriptional level, do not suggest the involvement of both transcription factors in the resistance to α-amanitin.

Genome Enrichment Analysis Confirms 30-Year-Old QTL Mapping Data

In order to identify the regulatory pathway components that lead to the α-amanitin resistance phenotype, we performed a genome enrichment analysis to look for clusters of significantly differentially expressed genes along the four chromosomes. In accordance with the two previous studies that mapped α-amanitin resistance to the polytene bands 95 and 98 on chromosome 3, we found signatures for both constitutive (group 2 versus group 1) and α-amanitin-inducible (group 3 versus group 2) clusters of differentially expressed genes. The only constitutively differentially expressed gene cluster is situated at cytological band 38B on the left arm of chromosome 2, which contains the genes CG10659, Taf13, CG17570, phr6-4, dia, and CG31674 at the peak of differential expression (Figure 3 and Table S4). However, their predicted and experimentally proven functions do not explain how α-amanitin resistance is genetically controlled. The remaining four clusters of differentially expressed genes responded to α-amanitin in the larval food. The most interesting induced gene cluster is situated at cytological band 66A, which is close to Mdr65-containing region 65A10 on the left arm of chromosome 3, to which α-amanitin resistance was previously mapped [14], [15]. The genes at the peak of differential expression are mp, Hsc70-4, pst, CG8562, Cyp316a1, Cyp4d8, CG33276, and RNaseX25 (Figure 3 and Table S4). Interestingly, two predicted Cytochrome P450 genes with unknown functions, Cyp316a1 and Cyp4d8, were 11.8- and 7.1-fold inducible by α-amanitin (see also Table 2). We further identified differentially expressed gene clusters at cytological bands 68A (left arm of chromosome 3), 92A, and 96D (right arm of chromosome 3). Most of these genes are poorly annotated and none of the genes were linked to any known toxin response (Figure 3 and Table S4). It is worth noting that the transcription factor gene Hr96 is close to the previously identified Pkc98 locus, to which α-amanitin resistance was mapped [14], [15]. Although our single gene analysis did not show significant up-regulation of the Hr96 gene, it is nevertheless possible that Hr96 contributes to the resistance on a post-transcriptional level.

[Figure omitted. See PDF.]

Figure 3. Genome enrichment analysis for genomic correlates.

Genomic correlates are likely disrupted in Ama-KTT/M/2 versus Canton S (red) and Ama-KTT/M/2 on α-amanitin versus Ama-KTT on non-toxic food (blue). Colored lines above the gray line indicate significant enrichment of a genomic correlate. Of the five genomic correlates rising above the cutoff value, two genomic correlates are similar to those found in previous linkage studies on the Ama-KTT stock.

https://doi.org/10.1371/journal.pone.0093489.g003

Gene Ontology Enrichment Analysis Suggests Additional α-Amanitin Resistance Mechanisms

In order to explore if multiple mechanisms confer the resistance phenotype to α-amanitin in the Ama-KTT/M/2 stock, we performed a gene ontology enrichment analysis. First, we compared the constitutive gene expression differences between Ama-KTT/M/2 and Canton-S on non-toxic food (group 2 versus group 1). As a result, we identified three molecular functions that could be relevant for the α-amanitin resistance phenotype in Ama-KTT/M/2 (Figure 4): 1) ‘Oxidoreductase activity’ genes (GO 0016491) were on average 4.6-fold higher expressed (p = 1.06E-18) in Ama-KTT/M/2. This result confirms the single gene analysis results (Table 1), which indicated that the three highest constitutively expressed Cyp genes (Cyp6a2, Cyp12d1-d, and Cyp12d1-p might be important for the resistance to α-amanitin. 2) ‘Transferase activity’ genes (GO 0016740) were on average 4.6-fold higher expressed in Ama-KTT/M/2 (p = 7.61E-11), confirming our single gene analysis for the Gst and Ugt genes (Table 1). 3) ‘Structural constituents of chitin-based cuticle’ genes (GO 0005214) were on average 10.5-fold (p = 1.87E-18) higher expressed in Ama-KTT/M/2, including 45 insect cuticle genes of the Cpr, Lcp, and Ccp gene families, which belong to the top 190 constitutively up-regulated genes in Ama-KTT/M/2 (Table S2). It is possible that cuticular proteins provide a protective layer against α-amanitin in organs that are covered by a cuticle, such as the epidermis and the gut. For example in honey bees, ‘structural constituents of chitin-based cuticle’ genes have been suggested to protect venom gland cells from toxins that are stored in the gland [49]. It is interesting to note that like α-amanitin, the bee venom ingredient Mast Cell Degranulating (MCD) Peptide is a bicyclic peptide. Structural constituents of the chitin-based cuticle could perhaps bind to bicyclic peptides and prevent them from entering cells. Furthermore, we identified two significant biological processes in this comparison (group 2 versus group 1) (Figure 4). 1) ‘Oxidation-reduction process’ genes (GO 0055114) were on average 5.6-fold higher expressed in Ama-KTT/M/2 (p = 5.25E-18), confirming the possible role of Cyp genes in α-amanitin detoxification. 2) The ‘cellular amino acid metabolic process’ genes (GO 0006520) showed a 1.2-fold higher expression average in Ama-KTT/M/2 (p = 2.55E-13) and was divided into two sub-processes. 2a) The ‘cellular modified amino acid process’ (GO 0006575) contained 16 Gst genes, which were on average 1.8-fold higher expressed in Ama-KTT/M/2 (p = 4.08E-03), suggesting that GST enzymes might help detoxifying α-amanitin via the phase II detoxification process. 2b) ‘Alpha-amino acid metabolic process’ genes (GO 1901605), such as glutathione metabolism genes, were on average 2.3-fold constitutively up-regulated in Ama-KTT/M/2 (p = 6.49E-03). Some of these genes might provide the substrate glutathione for the GST enzymes. Interestingly, yellow (y), a well-known pigmentation gene in Drosophila, was among the genes of this gene ontology term (14.7-fold up-regulated, p = 0.0448, Table S2). yellow is closely related to Major Royal Jelly Protein (MRJP) genes in honey bees, which were previously suggested to protect the venom gland cells from the bee venom [49]. It is thus possible that yellow plays a role in keeping α-amanitin outside of tissues or perhaps even modifying it so that it becomes less toxic.

[Figure omitted. See PDF.]

Figure 4. Gene ontology enrichment analysis for Ama-KTT/M/2 versus Canton-S on no toxin (group 2 versus 1).

The GO trees for the molecular function and biological process are shown on the left-hand side with the numbers for each term. The corrected p-values, average fold-changes for all genes in each term, term names, and selected genes of each GO term are shown on the right-hand side of each term number.

https://doi.org/10.1371/journal.pone.0093489.g004

Next, we aimed to identify the gene ontologies that respond to α-amanitin in the resistant stock Ama-KTT/M/2. We thus compared Ama-KTT/M/2 on α-amanitin-containing food to Ama-KTT/M/2 on non-toxic food (group 3 versus group 2). As a result, we identified genes with two molecular functions that are significantly induced by α-amanitin (Figure 5). 1) The ‘oxidoreductase activity’ genes (GO 0016491) are on average 4.7-fold induced (p = 2.36E-10) by the toxin, again suggesting that a phase I detoxification process mediated by Cytochrome P450s is involved in conferring α-amanitin resistance. Among the 37 Cyp genes of this gene ontology term, we found seven genes that we already identified in our single gene analysis (Table 2): Cyp316a1, Cyp6d2, Cyp4d8, Cyp28d1, Cyp6t1, Cyp4d2, and Cyp4d14. 2) ‘Peptidase activity, acting on L-amino acid peptides’ genes (GO 0070011) were on average 15.4-fold induced (p = 3.95E-05). Because α-amanitin is a peptide, peptidases are good candidates to cleave it. To date, however, no specific enzyme is known that can inactivate α-amanitin by cleaving this bicyclic octapeptide. Besides molecular functions, we further identified two biological processes that were of interest. 1) The ‘oxidation-reduction process’ genes (GO 0055114) were on average 5.0-fold induced (p = 3.40E-13), again confirming that Cyp genes could play a role in detoxifying α-amanitin. 2) We identified the ‘cellular amino acid metabolic process’ (GO 0006520) with an average up-regulation of 1.2-fold (p = 4.09E-11) in response to α-amanitin. The most interesting genes in this gene ontology group are 11 Gst genes and the yellow gene, again showing that the phase II detoxification process is inducible by α-amanitin and that yellow could play a protective role. Our gene ontology enrichment analysis further identified cellular components that respond to α-amanitin exposure (Figure 5). 1) ‘Cytoplasm’ genes (GO 0005737) were on average 681.2-fold induced (p = 5.26E-13), some of which are yellow, eight Cyp genes, and 13 Gst genes. The eight Cyp genes belong to the gene ontology term ‘cytoplasmic part’ (GO 0044444), which is on average 859.6-fold induced (p = 9.57E-10). Unexpectedly, the most highly induced gene ontology term for the cellular component was the ‘lipid particle’ with an average gene induction of 5,271.5-fold (p = 8.62E-10). Lipid particles are subcellular structures that play roles in detoxification processes and the innate immune system. In insects, lipid particles form coagulation products, thereby protecting cells from pathogens and toxic products of the phenol oxidase cascade [50]. In yeast cells, lipid particles detoxify excessive amounts of lipophilic substances [51]. Even in humans, liposomes are used for detoxifying patients with overdoses of drugs, such as heroin, opioids, and cocaine [52]. The fact that the Ama-KTT/M/2 stock responds to α-amanitin with a several thousand-fold induction of lipid particle genes suggests that cytoplasmic lipid particles contribute to the resistance to α-amanitin in the Ama-KTT/M/2 stock.

[Figure omitted. See PDF.]

Figure 5. Gene ontology enrichment analysis for Ama-KTT/M/2 on α-amanitin versus Ama-KTT/M/2 on no toxin (group 3 versus 2).

The GO trees for the molecular function, biological process, and cellular component are shown on the left-hand side with the numbers for each term. The corrected p-values, average fold-changes for all genes in each term, term names, and selected genes of each GO term are shown on the right-hand side of each term number.

https://doi.org/10.1371/journal.pone.0093489.g005

The Domain Enrichment Analysis Verifies the Gene Ontology Enrichment Analysis

Because many proteins have more than one functional domain and the gene ontology enrichment analysis cannot reveal what domain of a protein is important for the resistance to α-amanitin, we further performed a domain enrichment analysis with our microarray data. As shown in Table 3, when comparing Ama-KTT/M/2 with Canton-S (group 2 versus group 1) on non-toxic food, the following protein domains were identified as significantly enriched: Cytochrome P450 (p = 4.72E-11), UDP-glucuronosyl/UDP-glucosyltransferase (p = 1.26E-10), Cytochrome P450, conserved site (p = 5.93E-10), insect cuticle protein (p = 1.55E-09), Cytochrome P450, E-class, group I (p = 6.05E-09), Glutathione S-transferase, C-terminal (p = 6.39E-06), Glutathione S-transferase, C-terminal-like (p = 1.02E-05), Glutathione S-transferase/chloride channel, C-terminal (p = 1.28E-05), and Glutathione S-transferase, N-terminal (p = 4.77E-05). Thus, the domain enrichment analysis confirms the possible importance of phase I and II detoxification reactions in conferring α-amanitin resistance. When comparing Ama-KTT/M/2 on α-amanitin-containing food to Ama-KTT/M/2 on no toxin (group 3 versus group 2, Table 4), we identified the following significantly enriched protein domains: major royal jelly (p = 0), pupal cuticle protein C1 (p = 0), Cytochrome P450 (p = 1.20E-12), Cytochrome P450, conserved site (p = 2.90E-12), insect cuticle protein (p = 1.91E-11), chitin binding domain (p = 3.38E-11), Cytochrome P450, E-class, group I (p = 1.77E-10), peptidase M17, leucyl aminopeptidase, N-terminal (p = 4.38E-06), UDP-glucuronosyl/UDP-glucosyltransferase (p = 4.93E-06), leucine aminopeptidase/peptidase B (p = 8.48E-06), and peptidase M17, leucyl aminopeptidase, C-terminal (p = 8.49E-06). These results confirm the results from the gene ontology enrichment analysis, suggesting that Cytochrome P450s and transferases can detoxify α-amanitin via the phase I and II detoxification pathways. Furthermore, peptidases might cleave α-amanitin, and Royal Jelly Protein domain-containing proteins might protect tissues from α-amanitin, similar to the situation in the honey bee venom gland [49].

[Figure omitted. See PDF.]

Table 3. Domain enrichment analysis for Ama-KTT/M/2 versus Canton-S on no toxin (group 2 versus 1).

https://doi.org/10.1371/journal.pone.0093489.t003

[Figure omitted. See PDF.]

Table 4. Domain enrichment analysis for Ama-KTT/M/2 on α-amanitin versus Ama-KTT/M2 on no toxin (group 3 versus 2).

https://doi.org/10.1371/journal.pone.0093489.t004

The RT-qPCR Results Confirm the Microarray Data

We used real-time quantitative reverse transcription PCR (RT-qPCR) to confirm the fold-changes of ten genes, which we selected because of their high fold-changes and predicted importance for the resistance phenotype (Figure 6 and Table S5). When comparing Ama-KTT/M/2 to Canton-S (group 2 versus group 1), the genes Cyp6a2, 12d1-d, Ugt86Dd, GstD5, and GstE1 were between 1366.9 and 10.7-fold up-regulated (p<0.001 for all values, randomization test, B = 2000). When we compared Ama-KTT/M/2 treated with α-amanitin to Ama-KTT/M/2 (group 3 versus group 2), Cyp316a1, 6d2, 4d8, 28d1, and 6t1 were up-regulated between 14.1 and 8.4-fold (p = 0.002 for Cyp316a1 and p<0.001 for the other genes, randomization test, B = 2000). In summary, the microarray analysis fold-induction changes perfectly correlate with our RT-qPCR results, such that the microarray results slightly underestimate the fold-changes that resulted from the RT-qPCR analysis.

[Figure omitted. See PDF.]

Figure 6. The qPCR results confirm the microarray data.

A) Relative expression distribution (Y-axis) of ten selected genes is shown as a ratio comparing Ama-KTT/M/2 and Canton-S (group 2 versus group 1). Each measurement contains 15 replicates (3 replicates for each of the five biological controls of groups 1 and 2). B) Gene expression differences between Ama-KTT/M/2 treated with α-amanitin and Ama-KTT/M/2 (group 3 versus group 2) are compared. Group 3 contributes to 18 data points (three replicates for each of the six biological controls), while group 2 contributes to 15 data points, as previously mentioned. All comparisons were normalized with two reference genes, Sucb and alpha-Tub84B. Ratios above one indicate that a gene is up-regulated in the comparison.

https://doi.org/10.1371/journal.pone.0093489.g006

Discussion

Several Mechanisms Seem to Confer α-Amanitin Resistance

α-Amanitin is the principal toxin in some of the most deadly poisonous mushrooms, which inhibits the function of RNA-polymerase II by binding to it. Our results presented here comprise the first whole-transcriptome scale investigation to identify the molecular and cellular mechanisms that underlie the resistance to this very potent toxin in any organism. Using larvae of the resistant stock Ama-KTT/M/2 and the sensitive stock Canton-S, we identified both constitutive and α-amanitin-inducible mechanisms that can explain the resistance to α-amanitin in the Ama-KTT/M/2 stock. Based on an array of bioinformatics analyses of our microarray data and RT-qPCR validation, we found that four main mechanisms are likely to contribute in concert to the resistance: 1) constitutive and α-amanitin-inducible toxin entry blockage, mediated by cuticular proteins, the MRJP domain of the Yellow protein family, and Sgs proteins, 2) constitutive and α-amanitin-inducible phase I and II detoxification, mediated by the Cytochrome P450, GST, and UGT enzyme families (likely followed by excretion), 3) α-amanitin-inducible lipid particle gene induction, possibly leading to the sequestration of α-amanitin in cytoplasmic lipid particles, and 4) α-amanitin-inducible peptidase genes, perhaps leading to the digestion of α-amanitin either inside or outside (e.g. gut lumen) of cells (Figure 7).

[Figure omitted. See PDF.]

Figure 7. A model of the four mechanisms that contribute to the resistance to α-amanitin in concert.

The bicyclic octapeptide α-amanitin is shown as a red 8. Cuticular proteins block some of the α-amanitin from entering the cells (blockage). α-Amanitin that entered the cytoplasm is either sequestered in lipid particles, cleaved by peptidases, or detoxified by phase I and II detoxification enzymes, possibly followed by excretion.

https://doi.org/10.1371/journal.pone.0093489.g007

In honey bee venom glands, the Major Royal Jelly Protein 8 (MRJP8) was shown to be a part of the cuticular layer that forms the inner lining of the gland. It was suggested that MRJP8 protects the venom gland cells from the stored toxins [49]. The closest relatives to the MRJP genes in Drosophila are the proteins of the Yellow family. The yellow gene itself was together with numerous cuticular protein genes identified as significant in our single gene, gene ontology, and domain enrichment analyses. It is thus possible that Yellow, together with other cuticular proteins, block the entry of α-amanitin into cells protected by a cuticular layer, such as the larval epidermis and gut epithelium (Figure 7). In a similar manner, the products of the five strongly α-amanitin-inducible salivary gland secretion genes Sgs1, Sgs3, Sgs5, Sgs7, and Sgs8 (each >300-fold induced) could perhaps bind to α-amanitin and reduce its uptake in the midgut. Another possibility is that α-amanitin is simply a stress factor that induces these and other genes. After all, α-amanitin blocks messenger RNA transcription in poisoned cells, which is certainly stressful for the organism.

Besides being involved in environmental stress responses, hormone metabolism, and other metabolic functions, some Cytochrome P450, GST, and UGT proteins catalyze detoxification reactions, which transform a broad variety of xenobiotic substances into less toxic molecules that can be more easily excreted from the body [53]–[55]. Cytochrome P450 proteins, which are encoded by Cyp genes, are known for their broad range of substrates that they chemically modify. Several Cyp genes have been associated with single or multiple toxin resistance in diverse insect species, such as Cyp6g1 [5], [28], [56]–[58], Cyp6g2 [32], Cyp6a2 [18]–[25], Cyp12a4 [41], and Cyp12d1 [22], [23], [26]–[32]. Our single gene and gene ontology enrichment analyses identified three of these detoxification-implicated Cyp genes, which are more than about 200-fold constitutively up-regulated in Ama-KTT/M/2: Cyp6a2, Cyp12d1-d, and Cyp12d1-p (Table 1). It is thus possible that one or all three of these genes contribute to the resistance to α-amanitin. There is also evidence that Cyp12d1 is inducible by environmental stress factors, such as heat, oxidative stress, and air pollutants [29]–[31]. Because Cyp6a2, Cyp12d1-d, and Cyp12d1-p are constitutively up-regulated in our double-controlled study, stress is not a likely cause for the up-regulation of these three genes.

Some GST and UGT proteins perform phase II detoxification reactions that make toxic molecules bulkier and more hydrophobic, preparing the toxins for their excretion from the body. Several of these genes have been linked to insecticide resistance [6], [36], [54], [55], [59]–[71], while others are involved in several types of stress responses [34], [35], [45], [54]. Our single gene analysis showed that several Gst and Ugt genes are constitutively up-regulated in Ama-KTT-M/2 and that both gene families are significantly enriched in our gene ontology enrichment analysis, while their specific domains were identified as significant in the protein domain enrichment analysis. It is thus likely that some of them help detoxifying α-amanitin by making it both bulkier to prevent it from binding to RNA-Polymerase II and more water-soluble to augment its secretion via the Malpighian tubules (Figure 7). It is, however, possible that the α-amanitin-induced genes simply respond to stress caused by the effects of the toxin.

In our gene ontology enrichment analysis, we identified two other interesting mechanisms, which are inducible in response to α-amanitin in the larval food: the possible sequestration of α-amanitin in lipid particles and the cleavage of α-amanitin by peptidases. A group of genes involved in the cellular component ‘lipid particle’ were on average more than 5200-times induced by α-amanitin in the larval food. Natural and artificial lipid particles have been shown to be involved in various detoxification processes in very diverse organisms such as yeast, insects, and humans [50], [51], [72]. We therefore speculate that cytoplasmic lipid particles aggregate around α-amanitin molecules and trap them, thereby preventing the toxin from entering the nucleus, where RNA-Polymerase II performs its function. Furthermore, a variety of peptidase genes were identified in our various data analyses, suggesting that α-amanitin is cleaved either in the gut lumen, in the cells, or perhaps even in the food, if the larvae secrete peptidases from their mouths (Figure 7).

Implications

Our data does not support the previously held view that an MDR mechanism confers α-amanitin resistance in D. melanogaster. In 1982 and 2000, two studies based on QTL mapping suggested that α-amanitin-resistance in four wild-caught D. melanogaster stocks is conferred by two major loci on chromosome 3 [14], [15], the more recent of which pointed out Mdr65 and Pkc98E as possible candidates. However, our single gene and genome enrichment analyses identified two α-amanitin-inducible Cyp genes, Cyp316a1 and Cyp4d8, which are situated close to the Mdr65 locus and Hr96 close to the Pkc98E locus. Because Begun and Whitley used QTL mapping, not deletion mapping, the two Cyp and the Hr96 genes could instead be the resistance-conferring genes. Taking all the observations from our study together, we conclude that α-amanitin resistance has evolved as a quantitative complex trait that is based on entry blockage, phase I and II detoxification followed by secretion, peptidase cleavage, and sequestration.

Cross-resistance to a broad variety of toxins could explain how some Drosophila species evolved into mushroom-feeding specialists that can use mushroom toxins to their own advantage. For example, various mycophagous Drosophila species are frequently infected with parasitic nematodes that render about 20% of the adult flies sterile [11], [73]. Feeding on poisonous mushrooms not only kills the nematode parasites, it also provides a unique food source that is not accessible to many animals. D. melanogaster is a non-mycophagous species and should thus not be exposed to α-amanitin in nature. However, as discussed earlier, Cytochrome P450 enzymes can provide cross-resistance to multiple toxins, such as manufactured pesticides and natural xenobiotic products [32], [58]. We speculate that α-amanitin resistance in D. melanogaster has evolved in response to agricultural pesticides or other environmental factors, to which the flies were exposed before they were collected in the 1960s. Thus, if unrelated toxins can induce α-amanitin resistance, such a cross-resistance could prime a species to a radical host switch. If D. melanogaster females were to change their egg-laying behavior and oviposit on less toxic mushrooms, a niche change could result, followed by selection to feed on more toxic mushrooms. Being a species with such high fecundity, D. melanogaster could then even drive rare mycophagous Drosophila species out of their niche.

Limitations

The most obvious limitation of every microarray is that the observations and conclusions are entirely based on mRNA transcription differences. It is thus possible that some important mechanisms escaped detection. Furthermore, many D. melanogaster genes are still poorly annotated and their true functions are elusive. We thus excluded the most poorly annotated genes from our analysis. However, in doing so, we might have inadvertently lost some important genes that could contribute to the resistance to α-amanitin. Furthermore, because we used whole larvae in our study, we cannot determine the relative importance that the different tissues play in the resistance to α-amanitin.

Our microarray data analysis did not reveal any gene-regulatory pathways that lead to the resistance to α-amanitin. Hr96 and cnc have been shown to be upstream of detoxification genes [46], [47], [74]. Hr96 is situated on the right arm of the third chromosome, where the genome enrichment analysis shows a peak in response to α-amanitin. However, the expression levels of both Hr96 and cnc revealed nothing that would lead us to conclude their role in α-amanitin resistance. One reason for this could be that these genes encode transcription factors, which are already present in the cytoplasm to await activation, and we might not expect dramatic differences in their RNA regulation. Another reason could be that our larvae were feeding on α-amanitin from the first instar until they were collected at the late third instar. Thus, we might have missed the critical time period during which the upstream components of the pathway were up-regulated. We also noticed a lack of dramatic Cyp, Gst, and Ugt gene inducibility in response to α-amanitin. In the resistant stock Ama-KTT/M/2, many Cyp, Gst, and Ugt genes were constitutively expressed at higher levels than in Canton-S, while in larvae that were fed on toxic food, a completely different set of Cyp and Gst genes showed a much weaker induction than we initially expected. This weak gene induction is perhaps not surprising because in a previous microarray study using six different toxins, the detoxification gene families were not much inducible either [75]. It is thus possible that at least for the Cyp, Gst, and Ugt genes, the resistance to α-amanitin is mostly a constitutive trait.

Based on the mapping data from the two previous studies, we expected to find the α-amanitin resistance-conferring genes on chromosome 3 [14], [15]. Because the original Ama-KTT stock is 45 years old, we wanted to make sure that the genes on both major autosomes are homozygous before performing the microarray. One limitation to our approach is that we did not balance the X chromosome when we created the isochromosome stock Ama-KTT/M/2. However, we showed that the Ama-KTT/M/2 stock is not less resistant than original Ama-KTT stock (Figure 1), indicating that most if not all resistance-conferring alleles are present in the isochromosome stock that we used for the microarray. Most genes that we identified as significant are situated on chromosomes 2 and 3 (Tables 1 and 2). However, a few highly expressed genes, like yellow, are on the X chromosome. Thus, these X-chromosomal genes could either be the original alleles from Ama-KTT or the alleles from the multi-balancer stock. If they derived from the multi-balancer stock, the regulation of these genes could be explained by epistasis, such that the inducers of the X-chromosomal genes are situated on the two major autosomes, which are derived from the original Ama-KTT stock.

Future Studies

In order to identify the upstream components of the pathways that lead to the resistance to α-amanitin in the Ama-KTT/M/2 isochromosome stock, future microarray studies should include samples of larvae that have been exposed to α-amanitin for different periods of time. Because first instar larvae are very small, the exposure to α-amanitin should happen during the third larval instar, and samples should be collected at a series of subsequent time points thereafter. This approach should be efficient to detect gene-regulatory differences of the upstream pathway components. Furthermore, it would be interesting to investigate the mechanisms that confer α-amanitin resistance in mycophagous Drosophila species, using the RNA sequencing approach. Mycophagous species are several orders of magnitude more resistant to α-amanitin than D. melanogaster [11]–[13]. The higher toxin resistance of those species could produce clearer signals for the determination of the factors that make Drosophila resistant. After we gain a clearer picture about the candidate genes that might confer α-amanitin resistance in several Drosophila species, the next step would be to provide conclusive genetic evidence if the candidate genes are sufficient and necessary for the α-amanitin resistance phenotype. This could be done by transgenically overexpressing the resistance-conferring alleles in either D. melanogaster or other sensitive species that are closely related to highly resistant mycophagous species. In D. melanogaster, overexpression of candidate genes can be achieved using the Gal4-UAS system with visible effects in different organs such as the gut, fat body, and Malpighian tubules [32], [58]. Such tests can reveal the organs and tissues that contribute to the resistance to α-amanitin. Because toxic mushrooms contain more than one toxin, mycophagous Drosophila species must be resistant to a variety of toxins that target different biological processes [12], [76], [77]. Thus, other commercially available mushroom toxins, such as β-amanitin, phalloidin, ibotenic acid, and muscimol should be used to test if cross-resistance or independent mechanisms provide protection against the variety of mushroom toxins that mycophagous larvae encounter in their food source. Another pressing question is where α-amanitin goes once it entered a larva. Is it digested in the gut? Does it enter the cytoplasm of all or just a subset of cells? Radioactive α-amanitin could be a means to answer this question, but the analysis of the data might prove very difficult.

Conclusions

We suggest that the α-amanitin resistance phenotype in D. melanogaster, a species that does not feed on mushrooms in nature, has evolved as cross-resistance to pesticides or other factors in the environment. Entry blockage of α-amanitin into epithelial cells, phase I and II detoxification mediated by Cytochrome P450, GST, and UGT enzymes (likely to be followed by excretion from the body), sequestration of α-amanitin in cytoplasmic lipid particles, and proteolytic cleavage by peptidases are four likely mechanisms to contribute to the resistance phenotype in concert. In contrast, we did not detect any evidence for multidrug resistance efflux systems to be important for the resistance to α-amanitin. Future studies should include a time series of α-amanitin exposure, Drosophila species that actually feed on toxic mushrooms in nature, and more mushroom toxins. Candidate genes resulting from these experiments should then undergo sufficiency and necessity tests by transgenic rescue.

Materials and Methods

Fly Stocks

All fly stocks were maintained at room temperature on food containing Brewer's yeast, cornmeal, granulated sugar, agar, and methylparaben as antifungal agent. The wild-type stock Canton-S and the multi-balancer stock w[1118]/Dp(1;Y)y[+]; CyO/nub[1] b[1] sna[Sco] lt[1] stw[3]; MKRS/TM6B, Tb[1] were obtained from the Bloomington Stock Center, Bloomington, Indiana (stocks #1 and #3703, respectively). The α-amanitin-resistant Ama-KTT stock (# 14021-0231.07) was originally collected in 1968 in Kenting (Taiwan) and obtained from the Drosophila Species Stock Center at the University of California, San Diego.

Generation of the Isochromosome Stock Ama-KTT/M/2

Because Ama-KTT was maintained in the absence of selective pressure to toxins in the stock center over the past five decades, the stock could have lost, or become heterozygous for, some of the α-amanitin resistance-causing alleles. In order to create flies homozygous for the resistance-conferring alleles, we crossed the Ama-KTT stock to the multi-balancer stock w[1118]/Dp(1;Y)y[+]; CyO/nub[1] b[1] sna[Sco] lt[1] stw[3]; MKRS/TM6B, Tb[1]. As a result, we created the isochromosome stock Ama-KTT/M/2, which is isogenic for the second and third chromosomes.

Dose-Response Studies of the Fly Stocks to α-Amanitin

In order to quantify and compare the levels of α-amanitin resistance of the D. melanogaster stocks, dose-response experiments were performed, which measured the survival from freshly-hatched first-instar larvae to adulthood. Flies able to completely hatch from their pupae were scored as survivors. The α-amanitin-resistant stocks Ama-KTT and Ama-KTT/M/2 were tested on 11 α-amanitin concentrations, using 0 to 10 μg of α-amanitin per g of food in 1 μg increments. The α-amanitin-sensitive wild-type stocks Canton-S and the multi-balancer stock w[1118]/Dp(1;Y)y[+]; CyO/nub[1] b[1] sna[Sco] lt[1] stw[3]; MKRS/TM6B, Tb[1] were initially tested on five concentrations ranging from 0 to 4 μg of α-amanitin per g of food in 1 μg increments. However, because they survived only the zero-concentration, these stocks were further tested on 0, 0.25, 0.5, 0.75, 0.1, 0.25, and 0.375 μg of α-amanitin per g of food.

In order to obtain first-instar larvae for the dose-response experiments, flies of mixed sexes were allowed to lay eggs on molasses agar caps that contained a streak of fresh Baker's yeast paste at 25°C, 70% humidity, and a 12∶12 hour day/night cycle. The yeast was removed prior to larval hatching. Freshly hatched first-instar larvae were placed in groups of ten into 2-mL plastic test tubes (USA Scientific), each containing 500 mg of non-toxic or poisoned food and two small air holes in the lid. The food consisted of 125 mg dry, instant Drosophila medium (Carolina) and 375 μL sterile Milli-Q water with or without dissolved α-amanitin. Ten tubes were prepared for each toxin concentration and experimental replicate, resulting in 100 larvae for each concentration and experiment. Three high-quality dose-response experiments, in which the zero-concentration survival rate was at least 80%, were used to calculate the LC50 of each fly stock. The standard deviation of the mean (s.e.m.) was calculated for each concentration by sampling the data points of all 30 vials of every concentration. The LC50 was calculated using scatter plots and the logarithmic trendline function in Microsoft Excel.

Sample Preparation for the Microarray Analysis

In order to compare the constitutive gene-regulatory differences across the entire transcriptome between α-amanitin-sensitive and -resistant stocks, freshly-hatched first-instar larvae of the sensitive Canton-S stock (group 1) and the resistant Ama-KTT/M/2 stock (group 2) were placed in groups of ten into 2-mL plastic test tubes (USA Scientific), containing 500 mg of non-toxic food. To identify the genes that are inducible by α-amanitin, Ama-KTT/M/2 larvae were raised on 1.5 μg of α-amanitin per g of food (group 3), which is slightly lower than the LC50 concentration of this stock. All larvae were raised until they reached the late third instar at 25°C, 70% humidity, and a 12∶12 hour day/night cycle. Because not all larvae survived in the tubes and the larvae on α-amanitin-containing food had a slower growth rate, initially 600 first-instar larvae (60 tubes) for each group were started over three subsequent days (20 tubes per group and day). When the majority of larvae reached the late third instar, the tubes were emptied and groups of ten late, but still feeding third-instar larvae were randomly picked from across all tubes and flash-frozen in batches of ten in liquid nitrogen, each batch providing the RNA for one microarray chip. Five biological replicates (ten larvae each) were prepared for groups 1 and 2, whereas group 3 was prepared in six biological replicates (ten larvae each). All samples were collected on the same morning. RNA extraction was performed without delay, using the RNeasy microarray tissue kit (Qiagen), according to the manufacturer's instructions.

Affymetrix Array Target Preparation, Hybridization, and Scanning

Collection and analysis of data were compliant with MIAME standards [78]. The microarray experiment was performed using the Affymetrix GeneChip Drosophila Genome 2.0 Arrays (Affymetrix, Santa Clara, CA, USA) with biotinylated targets derived from total RNA. Each array contains 18,952 probes that interrogate ∼18500 transcripts of genes present in the transcriptome of D. melanogaster. Prior to labeling, total RNA samples were checked for purity and concentration, using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher, Waltham, MA, USA) and for integrity, using RNA 6000 Nano Chips in a BioAnalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). cDNA for hybridization was synthesized and biotin-labeled from 400 ng of total RNA, using a MessageAmp Premier IVT kit (Ambion, Austin, TX, USA) according to the manufacturer's specifications. Biotinylated cDNA was fragmented, then hybridized, washed, and stained using a GeneChip Hybridization, Wash, and Stain Kit (Affymetrix, Santa Clara, CA, USA) according to the manufacturer's specifications. Arrays were post-processed on the AFX 450 Fluidics Station before they were scanned on an AFX GC3000 G7 Scanner (Affymetrix, Austin, TX, USA). Data were extracted from the raw images, using the Affymetrix Expression Console v.1.2 software. The RNA quality check, labeling, hybridization, and imaging procedures were performed according to Affymetrix protocols at the Center for Genomics Research and Biocomputing, University of Wisconsin.

Microarray Data Normalization

The quality of microarray data sets was first checked by examining the distribution of the Studentized deleted residuals, using a previously described procedure [79], [80], and only high-quality microarray data were used for normalization. Probeset-level normalization was performed with the PLIER (Probe Logarithmic Intensity Error) algorithm with quantile normalization and mismatch intensity adjustment, using the Affymetrix Power Tools software v.1.14.4.1. Probesets were annotated using release 32 of the Affymetrix annotation for the Drosophila 2.0 array platform. The CEL files and summarized (normalized) microarray data resulting from this study have been deposited in the NCBI's Gene Expression Omnibus database at NIH (http://www.ncbi.nlm.nih.gov/geo/) with the accession number of GSE52782.

Genome Enrichment Analysis

To find genome regions containing more differentially expressed genes than expected by chance, we used the binomial coincidence detection algorithm [81] with modifications specific for this dataset. Because D. melanogaster has a smaller genome and shorter regions of genetic linkage than mammalian genomes, we reduced the length of the overlapping bins to 500 kb spaced at 250 kb intervals. In order to reduce the total noise and find the strongest signal, we used only the top 0.01 most differentially expressed genes in the dataset. Briefly, under a null hypothesis of no significant enrichment in a genome region, the probability of finding a significantly differentially expressed gene within each bin will follow a binomial distribution with a probability of any given gene being significantly differentially expressed at no more than 0.01. The algorithm calculates a binomial probability for the empirical quantity of differentially expressed genes within each bin across the entire genome. The decimal log of the inverse of these probabilities is graphed. A decimal log of 2 corresponding to the horizontal line through each graph indicates a probability of a cluster occurring 1 in 100 times under the null hypothesis, the cutoff used for this method. The resulting graphic shows clustering over the whole genome and spikes indicate clusters unlikely to have occurred by chance. This is statistical evidence that a genome region is likely implicated in a gene expression phenotype. The assumptions for the inferential statistics used for this analysis necessitate inclusion of low copy genes as differentially expressed, thus the inferential statistics used to generate the genome enrichment figure were performed in the limma package in Bioconductor v.2.10 [82]. Cytoband visualization is derived from annotation tables of the UCSC dm3 genome, which represents cytobands as alternating light and dark bands.

Identification of Differentially Expressed Genes (DEGs)

A nonparametric method, RankProd (RP) [83], was used to identify differentially expressed genes (DEGs) between Ama-KTT treated with α-amanitin slightly below the LC50 concentration, and untreated Ama-KTT, or Canton-S conditions. We chose RP because it had been implicated to be more accurate for ranking genes by differential expression than t-statistics or derived methods [84]. Kadota et al. once evaluated eight DEG ranking methods and concluded that RP is one of the best performing methods [84]. Laing et al. indicated RP is one of most efficient method when replicate numbers is less than 10 [85]. In this study, we applied multiple testing corrections to the p-values resulting from RP using Benjamini and Hochberg False Discovery Rate [86] and all genes with corrected p-values (<0.05) were defined as DEGs.

Gene Ontology Enrichment Analysis

The DEGs identified from each comparison, namely, Ama-KTT/M/2 versus Canton-S and Ama-KTT/M/2 on α-amanitin versus Ama-KTT/M/2 on no toxin, were used as the input for the gene ontology enrichment analysis. We employed an online tool, AmiGO's Term Enrichment, to identify the enriched gene ontologies (http://amigo.geneontology.org/). This tool uses the Perl module GO:TermFinder available at CPAN (http://search.cpan.org/) to identify the enriched gene ontology terms associated with a DEG list, using the hypergeometric probability function. We applied multiple testing corrections to calculate the p-values of all GO terms and then corrected p-values using Benjamini and Hochberg False Discovery Rate [86]. All gene ontology terms with a corrected p-value <0.05 were considered to be significantly enriched.

Protein Domain Enrichment (PDE) Analysis

Protein domains were analyzed with InterproScan [87]. We first downloaded and installed InterproScan and associated databases to our Linux server and performed the standalone analysis to identify protein domains of all target sequences provided by FlyBase (http://flybase.org/static_pages/docs/datafiles.html). The enrichment of each domain in the differentially expressed gene list was compared to the occurrence of the respective domain in the background of all genomic genes, and two parameters were introduced to show the enrichment of each domain as described in [88]: (1) Enrichment factor, EF = k/(nM/N); and (2) the E_score, which is the hypergeometric probability of identifying at least k domains from DEG list. It is calculated using the following formula:

N is the total number of domains associated with all genomic genes, M is total number of a specific domain for all genes in the genome, n is the number of all domains associated with the DEGs, and k is the number of a specific domain present in the DEGs list. We applied multiple testing corrections to the p-values calculated via hypergeometric probability using Benjamini and Hochberg False Discovery Rate (FDR) [86]. The significantly enriched protein domains are those that have a corrected p-value <0.05.

RT-qPCR Validation of the Microarray Results

Quantitative real-time PCR (qPCR) was performed on ten genes of interest to confirm the results of the microarray analysis. Each gene of interest and biological replicate was repeated three times to ensure the statistical significance of the result. The genes included two Cyp, one Ugt, and two Gst genes that were up-regulated when comparing the resistant group Ama-KTT/M/2 to the control group Canton-S (group 2 versus group 1) and five Cyp genes that were up-regulated when comparing Ama-KTT/M/2 on α-amanitin to Ama-KTT/M/2 on no toxin (group 3 versus group 2). Two reference genes, Scub and alpha-tub84B, were used as controls to normalize the results. These genes were selected because their fold-changes were nearly zero for each comparison. The primer pairs used were a part of the Taqman Gene Expression Assays kit (Applied Biosystems): Dm02361072_s1, Dm01831596_g1, Dm01840671_g1, Dm01830394_g1, Dm01822311_g1, Dm01804633_g1, Dm01799869_s1, Dm02147253_g1, Dm01817955_g1, Dm02152265_s1, Dm01826948_s1, and Dm02374415_g1. The reactions were performed in a StepOnePlus Real-Time PCR System (Applied Biosystems). The High Capacity cDNA Reverse Transcription Kit (Applied Biosystems) was used to reverse transcribe RNA to cDNA in an Eppendorf PCR machine for 96 reactions (Eppendorf, Model 96S). We used REST 2009 to calculate the RT-qPCR p-values.

Supporting Information

[Figure omitted. See PDF.]

Table S1.

Differentially expressed genes (DEGs) between Ama-KTT/M/2 on no toxin versus Canton-S (group 2 versus 1), Ama-KTT/M/2 on toxin versus Canton-S (group 3 versus 1), and Ama-KTT/M/2 on toxin versus Ama-KTT/M/2 on no toxin (group 3 versus 2). This table contains 4209 DEGs that are differentially expressed in at least one of the three comparisons.

https://doi.org/10.1371/journal.pone.0093489.s001

(XLSX)

Table S2.

Complete single gene analysis for Ama-KTT/M/2 versus Canton-S on no toxin (group 2 versus 1). This table contains well-annotated genes that are at least 2.0-fold constitutively up-regulated in the resistant stock, as compared to the sensitive stock, on no toxin. The p-value cutoff is p<0.05.

https://doi.org/10.1371/journal.pone.0093489.s002

(XLSX)

Table S3.

Complete single gene analysis for Ama-KTT/M/2 on toxin versus Ama-KTT/M/2 on no toxin (group 3 versus 2). This table contains well-annotated genes that are at least 2.0-fold inducible by feeding larvae of the resistant stock with α-amanitin-containing food, as compared to resistant larvae on no toxin. The p-value cutoff is p<0.05.

https://doi.org/10.1371/journal.pone.0093489.s003

(XLSX)

Table S4.

Genome enrichment analysis for group 2 versus 1 and group 3 versus 2. This table shows the genes behind the peaks in Figure 3. The peak at band 38B is the only locus that is differentially expressed between Ama-KTT/M/2 and Canton-S on no toxin (group 2 versus group 1). The remaining peaks 66A, 69A, 92A, and 96D show differentially expressed loci in response to α-amanitin treatment (group 3 versus group 2). All p-values are corrected and fold-changes are given for the individual genes. Peaks 66A and 96D are very close to the two QTL mapping peaks identified in previous studies [14], [15].

https://doi.org/10.1371/journal.pone.0093489.s004

(XLSX)

Table S5.

Comparison of qPCR and microarray fold-induction values. The first five genes were constitutively over-expressed in Ama-KTT/M/2, as Compared to Canton-S (Group 2 versus group 1). The last five genes were induced by α-amanitin in Ama-KTT/M/2, as compared to Ama-KTT/M/2 on no toxin (group 3 versus group 2). The RT-qPCR p-values are uncorrected, while the array p-values are corrected.

https://doi.org/10.1371/journal.pone.0093489.s005

(XLSX)

Acknowledgments

We thank Felicia Nip, Zachary Johnson, and Kara Vogel for technical assistance with the sample collection. We would also like to thank Wayne Davis at the University of Wisconsin-Madison Biotechnology Center for his work on the microarray project and John Jaenike at the University of Rochester, NY, for helpful discussions during the preparation of the manuscript.

Author Contributions

Conceived and designed the experiments: CM TW. Performed the experiments: CM TW. Analyzed the data: CM MS LL HW TW. Contributed reagents/materials/analysis tools: MS HW TW. Wrote the paper: CM MS LL HW TW.

Citation: Mitchell CL, Saul MC, Lei L, Wei H, Werner T (2014) The Mechanisms Underlying α-Amanitin Resistance in Drosophila melanogaster: A Microarray Analysis. PLoS ONE 9(4): e93489. https://doi.org/10.1371/journal.pone.0093489

References

1. Despres L, David JP, Gallet C (2007) The evolutionary ecology of insect resistance to plant chemicals. Trends Ecol Evol 22: 298–307.

2. Feldman CR, Brodie ED Jr, Brodie ED 3rd, Pfrender ME (2012) Constraint shapes convergence in tetrodotoxin-resistant sodium channels of snakes. Proc Natl Acad Sci U S A 109: 4556–4561.

3. Bricelj VM, Connell L, Konoki K, Macquarrie SP, Scheuer T, et al. (2005) Sodium channel mutation leading to saxitoxin resistance in clams increases risk of PSP. Nature 434: 763–767.

4. Baxter SW, Badenes-Perez FR, Morrison A, Vogel H, Crickmore N, et al. (2011) Parallel evolution of Bacillus thuringiensis toxin resistance in lepidoptera. Genetics 189: 675–679.

5. Daborn PJ, Yen JL, Bogwitz MR, Le Goff G, Feil E, et al. (2002) A single P450 allele associated with insecticide resistance in Drosophila. Science 297: 2253–2256.

6. Enayati AA, Ranson H, Hemingway J (2005) Insect glutathione transferases and insecticide resistance. Insect Mol Biol 14: 3–8.

7. Ffrench-Constant RH (2013) The molecular genetics of insecticide resistance. Genetics 194: 807–815.

8. Ranson H, Claudianos C, Ortelli F, Abgrall C, Hemingway J, et al. (2002) Evolution of supergene families associated with insecticide resistance. Science 298: 179–181.

9. Pimentel D, Acquay H, Biltonen M, Rice P, Silva M, et al. (1992) Environmental and Economic Costs of Pesticide Use. Bioscience 42: 750–760.

10. Lindell TJ, Weinberg F, Morris PW, Roeder RG, Rutter WJ (1970) Specific inhibition of nuclear RNA polymerase II by alpha-amanitin. Science 170: 447–449.

11. Jaenike J (1985) Parasite Pressure and the Evolution of Amanitin Tolerance in Drosophila. Evolution 39: 1295–1301.

12. Stump AD, Jablonski SE, Bouton L, Wilder JA (2011) Distribution and Mechanism of alpha-Amanitin Tolerance in Mycophagous Drosophila (Diptera: Drosophilidae). Environ Entomol 40: 1604–1612.

13. Jaenike J, Grimaldi DA, Sluder AE, Greenleaf AL (1983) Alpha-Amanitin Tolerance in Mycophagous Drosophila. Science 221: 165–167.

14. Phillips JP, Willms J, Pitt A (1982) alpha-amanitin resistance in three wild strains of Drosophila melanogaster. Can J Genet Cytol 24: 151–162.

15. Begun DJ, Whitley P (2000) Genetics of alpha-amanitin resistance in a natural population of Drosophila melanogaster. Heredity 85: 184–190.

16. Chambers TC, McAvoy EM, Jacobs JW, Eilon G (1990) Protein kinase C phosphorylates P-glycoprotein in multidrug resistant human KB carcinoma-cells. J Biol Chem 265: 7679–7686.

17. Chung H, Sztal T, Pasricha S, Sridhar M, Batterham P, et al. (2009) Characterization of Drosophila melanogaster cytochrome P450 genes. Proc Natl Acad Sci U S A 106: 5731–5736.

18. Brun A, Cuany A, Le Mouel T, Berge J, Amichot M (1996) Inducibility of the Drosophila melanogaster cytochrome P450 gene, CYP6A2, by phenobarbital in insecticide susceptible or resistant strains. Insect Biochem Mol Biol 26: 697–703.

19. Dunkov BC, Guzov VM, Mocelin G, Shotkoski F, Brun A, et al. (1997) The Drosophila cytochrome P450 gene Cyp6a2: structure, localization, heterologous expression, and induction by phenobarbital. DNA Cell Biol 16: 1345–1356.

20. Giraudo M, Unnithan GC, Le Goff G, Feyereisen R (2010) Regulation of cytochrome P450 expression in Drosophila: Genomic insights. Pestic Biochem Physiol 97: 115–122.

21. Kalajdzic P, Oehler S, Reczko M, Pavlidi N, Vontas J, et al. (2012) Use of mutagenesis, genetic mapping and next generation transcriptomics to investigate insecticide resistance mechanisms. PLoS One 7: e40296.

22. Pedra JHF, McIntyre LM, Scharf ME, Pittendrigh BR (2004) Genom-wide transcription profile of field- and laboratory-selected dichlorodiphenyltrichloroethane (DDT)-resistant Drosophila. Proc Natl Acad Sci U S A 101: 7034–7039.

23. Sun W, Margam VM, Sun L, Buczkowski G, Bennett GW, et al. (2006) Genome-wide analysis of phenobarbital-inducible genes in Drosophila melanogaster. Insect Mol Biol 15: 455–464.

24. Wan H, Liu Y, Li M, Zhu S, Li X, et al. (2013) Nrf2/Maf binding site-containing functional Cyp6a2 allele is associated with DDT resistance in Drosophila melanogaster. Pest Manag Sci.

25. Amichot M, Tares S, Brun-Barale A, Arthaud L, Bride JM, et al. (2004) Point mutations associated with insecticide resistance in the Drosophila cytochrome P450 Cyp6a2 enable DDT metabolism. Eur J Biochem 271: 1250–1257.

26. Brandt A, Scharf M, Pedra JH, Holmes G, Dean A, et al. (2002) Differential expression and induction of two Drosophila cytochrome P450 genes near the Rst(2)DDT locus. Insect Mol Biol 11: 337–341.

27. Festucci-Buselli RA, Carvalho-Dias AS, de Oliveira-Andrade M, Caixeta-Nunes C, Li HM, et al. (2005) Expression of Cyp6g1 and Cyp12d1 in DDT resistant and susceptible strains of Drosophila melanogaster. Insect Mol Biol 14: 69–77.

28. Le Goff G, Boundy S, Daborn PJ, Yen JL, Sofer L, et al. (2003) Microarray analysis of cytochrome P450 mediated insecticide resistance in Drosophila. Insect Biochem Mol Biol 33: 701–708.

29. McDonnell CM, King D, Comeron JM, Li H, Sun W, et al. (2012) Evolutionary toxicogenomics: diversification of the Cyp12d1 and Cyp12d3 genes in Drosophila species. J Mol Evol 74: 281–296.

30. Moskalev A, Shaposhnikov M, Snezhkina A, Kogan V, Plyusnina E, et al. (2014) Mining gene expression data for pollutants (dioxin, toluene, formaldehyde) and low dose of gamma-irradiation. PLoS One 9: e86051.

31. Sun LJ, Schemerhorn B, Jannasch A, Walters KR, Adamec J, et al. (2011) Differential transcription of cytochrome P450s and glutathione S transferases in DDT-susceptible and -resistant Drosophila melanogaster strains in response to DDT and oxidative stress. Pestic Biochem Physiol 100: 7–15.

32. Daborn PJ, Lumb C, Boey A, Wong W, Ffrench-Constant RH, et al. (2007) Evaluating the insecticide resistance potential of eight Drosophila melanogaster cytochrome P450 genes by transgenic over-expression. Insect Biochem Mol Biol 37: 512–519.

33. Guzov VM, Unnithan GC, Chernogolov AA, Feyereisen R (1998) CYP12A1, a mitochondrial cytochrome P450 from the house fly. Arch Biochem Biophys 359: 231–240.

34. Yepiskoposyan H, Egli D, Fergestad T, Selvaraj A, Treiber C, et al. (2006) Transcriptome response to heavy metal stress in Drosophila reveals a new zinc transporter that confers resistance to zinc. Nucleic Acids Res 34: 4866–4877.

35. Neal SJ, Karunanithi S, Best A, So AKC, Tanguay RM, et al. (2006) Thermoprotection of synaptic transmission in a Drosophila heat shock factor mutant is accompanied by increased expression of Hsp83 and DnaJ-1. Physiol Genomics 25: 493–501.

36. Lumjuan N, Rajatileka S, Changsom D, Wicheer J, Leelapat P, et al. (2011) The role of the Aedes aegypti Epsilon glutathione transferases in conferring resistance to DDT and pyrethroid insecticides. Insect Biochem Mol Biol 41: 203–209.

37. Dierick HA, Greenspan RJ (2006) Molecular analysis of flies selected for aggressive behavior. Nat Genet 38: 1023–1031.

38. Huang QY, Sun PD, Zhou XG, Lei CL (2012) Characterization of Head Transcriptome and Analysis of Gene Expression Involved in Caste Differentiation and Aggression in Odontotermes formosanus (Shiraki). Plos One 7: e50383

39. Wang LM, Dankert H, Perona P, Anderson DJ (2008) A common genetic target for environmental and heritable influences on aggressiveness in Drosophila. Proc Natl Acad Sci U S A 105: 5657–5663.

40. Hansen BH, Altin D, Hessen KM, Dahl U, Breitholtz M, et al. (2008) Expression of ecdysteroids and cytochrome P450 enzymes during lipid turnover and reproduction in Calanus finmarchicus (Crustacea: Copepoda). General and Comparative Endocrinology 158: 115–121.

41. Bogwitz MR, Chung H, Magoc L, Rigby S, Wong W, et al. (2005) Cyp12a4 confers lufenuron resistance in a natural population of Drosophila melanogaster. Proc Natl Acad Sci U S A 102: 12807–12812.

42. Wang SP, He GL, Chen RR, Li F, Li GQ (2012) The Involvement of Cytochrome P450 Monooxygenases in Methanol Elimination in Drosophila Melanogaster Larvae. Arch Insect Biochem Physiol 79: 264–275.

43. Kang J, Kim J, Choi KW (2011) Novel cytochrome P450, cyp6a17, is required for temperature preference behavior in Drosophila. PLoS One 6: e29800.

44. Thomas AM, Hui C, South A, McVey M (2013) Common Variants of Drosophila melanogaster Cyp6d2 Cause Camptothecin Sensitivity and Synergize With Loss of Brca2. G3 3: 91–99.

45. Jumbo-Lucioni PP, Hopson ML, Hang D, Liang YL, Jones DP, et al. (2013) Oxidative stress contributes to outcome severity in a Drosophila melanogaster model of classic galactosemia. Dis Model Mech 6: 84–94.

46. King-Jones K, Horner MA, Lam G, Thummel CS (2006) The DHR96 nuclear receptor regulates xenobiotic responses in Drosophila. Cell Metab 4: 37–48.

47. Sykiotis GP, Bohmann D (2008) Keap1/Nrf2 signaling regulates oxidative stress tolerance and lifespan in Drosophila. Dev Cell 14: 76–85.

48. Sykiotis GP, Bohmann D (2010) Stress-Activated Cap'n'collar Transcription Factors in Aging and Human Disease. Sci Signal 3.

49. Peiren N, de Graaf DC, Vanrobaeys F, Danneels EL, Devreese B, et al. (2008) Proteomic analysis of the honey bee worker venom gland focusing on the mechanisms of protection against tissue damage. Toxicon 52: 72–83.

50. Rahman MM, Ma G, Roberts HLS, Schmidt O (2006) Cell-free immune reactions in insects. J Insect Physiol 52: 754–762.

51. Mullner H, Daum G (2004) Dynamics of neutral lipid storage in yeast. Acta Biochim Pol 51: 323–347.

52. Howell BA, Chauhan A (2010) Current and Emerging Detoxification Therapies for Critical Care. Materials 3: 2483–2505.

53. Mansuy D (1998) The great diversity of reactions catalyzed by cytochromes P450. Comp Biochem Physiol C Pharmacol Toxicol Endocrinol 121: 5–14.

54. Ranson H, Hemingway J (2005) Mosquito glutathione transferases. Gluthione Transferases and Gamma-Glutamyl Transpeptidases 401: 226–241.

55. Luque T, Okano K, O'Reilly DR (2002) Characterization of a novel silkworm (Bombyx mori) phenol UDP-glucosyltransferase. Eur J Biochem 269: 819–825.

56. Daborn P, Boundy S, Yen J, Pittendrigh B, ffrench-Constant R (2001) DDT resistance in Drosophila correlates with Cyp6g1 over-expression and confers cross-resistance to the neonicotinoid imidacloprid. Mol Genet Genomics 266: 556–563.

57. Joussen N, Heckel DG, Haas M, Schuphan I, Schmidt B (2008) Metabolism of imidacloprid and DDT by P450 CYP6G1 expressed in cell cultures of Nicotiana tabacum suggests detoxification of these insecticides in Cyp6g1-overexpressing strains of Drosophila melanogaster, leading to resistance. Pest Manag Sci 64: 65–73.

58. Harrop TW, Sztal T, Lumb C, Good RT, Daborn PJ, et al. (2014) Evolutionary changes in gene expression, coding sequence and copy-number at the cyp6g1 locus contribute to resistance to multiple insecticides in Drosophila. PLoS One 9: e84879.

59. Clark AG, Shamaan NA (1984) Evidence That Ddt-Dehydrochlorinase from the Housefly Is a Glutathione S-Transferase. Pestic Biochem Physiol 22: 249–261.

60. Clark AG, Shamaan NA, Sinclair MD, Dauterman WC (1986) Insecticide Metabolism by Multiple Glutathione S-Transferases in 2 Strains of the Housefly, Musca-Domestica (L). Pestic Biochem Physiol 25: 169–175.

61. Wang JY, Mccommas S, Syvanen M (1991) Molecular-Cloning of a Glutathione-S-Transferase Overproduced in an Insecticide-Resistant Strain of the Housefly (Musca-Domestica). Mol Gen Genet 227: 260–266.

62. Fournier D, Bride JM, Poirie M, Berge JB, Plapp FW (1992) Insect Glutathione S-Transferases - Biochemical Characteristics of the Major Forms from Houseflies Susceptible and Resistant to Insecticides. J Biol Chem 267: 1840–1845.

63. Prapanthadara LA, Hemingway J, Ketterman AJ (1993) Partial-Purification and Characterization of Glutathione S-Transferases Involved in Ddt Resistance from the Mosquito Anopheles-Gambiae. Pestic Biochem Physiol 47: 119–133.

64. Prapanthadara LA, Koottathep S, Promtet N, Hemingway J, Ketterman AJ (1996) Purification and characterization of a major glutathione S-transferase from the mosquito Anopheles dirus (Species b). Insect Biochem Mol Biol 26: 277–285.

65. Penilla RP, Rodriguez AD, Hemingway J, Torres JL, Arredondo-Jimenez JI, et al. (1998) Resistance management strategies in malaria vector mosquito control. Baseline data for a large-scale field trial against Anopheles albimanus in Mexico. Med Vet Entomol 12: 217–233.

66. Wei SH, Clark AG, Syvanen M (2001) Identification and cloning of a key insecticide-metabolizing glutathione S-transferase (MdGST-6A) from a hyper insecticide-resistant strain of the housefly Musca domestica. Insect Biochem Mol Biol 31: 1145–1153.

67. Ranson H, Rossiter L, Ortelli F, Jensen B, Wang XL, et al. (2001) Identification of a novel class of insect glutathione S-transferases involved in resistance to DDT in the malaria vector Anopheles gambiae. Biochem J 359: 295–304.

68. Ortelli F, Rossiter LC, Vontas J, Ranson H, Hemingway J (2003) Heterologous expression of four glutathione transferase genes genetically linked to a major insecticide-resistance locus from the malaria vector Anopheles gambiae. Biochem J 373: 957–963.

69. Lumjuan N, McCarroll L, Prapanthadara LA, Hemingway J, Ranson H (2005) Elevated activity of an Epsilon class glutathione transferase confers DDT resistance in the dengue vector, Aedes aegypti. Insect Biochem Mol Biol 35: 861–871.

70. Alias Z, Clark AG (2007) Studies on the glutathione S-transferase proteome of adult Drosophila melanogaster: responsiveness to chemical challenge. Proteomics 7: 3618–3628.

71. Gunasekaran K, Muthukumaravel S, Sahu SS, Vijayakumar T, Jambulingam P (2011) Glutathione S transferase activity in Indian vectors of malaria: A defense mechanism against DDT. J Med Entomol 48: 561–569.

72. Howell B, Chauhan A (2008) Uptake of amitriptyline and nortriptyline with liposomes, proteins, and serum: implications for drug detoxification. J Colloid Interface Sci 319: 81–93.

73. Perlman SJ, Jaenike J (2003) Infection success in novel hosts: an experimental and phylogenetic study of Drosophila-parasitic nematodes. Evolution 57: 544–557.

74. Misra JR, Horner MA, Lam G, Thummel CS (2011) Transcriptional regulation of xenobiotic detoxification in Drosophila. Genes Dev 25: 1796–1806.

75. Willoughby L, Chung H, Lumb C, Robin C, Batterham P, et al. (2006) A comparison of Drosophila melanogaster detoxification gene induction responses for six insecticides, caffeine and phenobarbital. Insect Biochem Mol Biol 36: 934–942.

76. Hallen HE, Luo H, Scott-Craig JS, Walton JD (2007) Gene family encoding the major toxins of lethal Amanita mushrooms. Proc Natl Acad Sci U S A 104: 19097–19101.

77. Tuno N, Takahashi KH, Yamashita H, Osawa N, Tanaka C (2007) Tolerance of Drosophila flies to ibotenic acid poisons in mushrooms. J Chem Ecol 33: 311–317.

78. Brazma A, Hingamp P, Quackenbush J, Sherlock G, Spellman P, et al. (2001) Minimum information about a microarray experiment (MIAME) - toward standards for microarray data. Nat Genet 29: 365–371.

79. Persson S, Wei HR, Milne J, Page GP, Somerville CR (2005) Identification of genes required for cellulose synthesis by regression analysis of public microarray data sets. Proc Natl Acad Sci U S A 102: 8633–8638.

80. Wei H, Persson S, Mehta T, Srinivasasainagendra V, Chen L, et al. (2006) Transcriptional coordination of the metabolic network in Arabidopsis. Plant Physiol 142: 762–774.

81. Saul MC, Gessay GM, Gammie SC (2012) A new mouse model for mania shares genetic correlates with human bipolar disorder. PLoS One 7: e38128.

82. Smyth GK (2004) Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol 3: Article3.

83. Breitling R, Armengaud P, Amtmann A, Herzyk P (2004) Rank products: a simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments. FEBS Lett 573: 83–92.

84. Kadota K, Shimizu K (2011) Evaluating methods for ranking differentially expressed genes applied to microArray quality control data. BMC Bioinformatics 12: 227.

85. Laing E, Smith CP (2010) RankProdIt: A web-interactive Rank Products analysis tool. BMC Res Notes: 221.

86. Benjamini Y, Hochberg Y (1995) Controlling the False Discovery Rate - a Practical and Powerful Approach to Multiple Testing. J R Stat Soc Series B Stat Methodol 57: 289–300.

87. Zdobnov EM, Apweiler R (2001) InterProScan - an integration platform for the signature-recognition methods in InterPro. Bioinformatics 17: 847–848.

88. Wei HR, Gou JQ, Yordanov Y, Zhang HX, Thakur R, et al. (2013) Global transcriptomic profiling of aspen trees under elevated [CO2] to identify potential molecular mechanisms responsible for enhanced radial growth. J Plant Res 126: 305–320.