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

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.


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][5][6][7][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 a-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 aamanitin-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 a-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 a-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 aamanitin than the sensitive wild-type strain Oregon-R [14]. In these three Asian strains, the resistance to a-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 aamanitin 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 a-amanitin from cells, the question of how D. melanogaster evolved a-amanitin resistance appeared to be answered. Although the proposed scenario is simple and elegant, no conclusive evidence has been

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 a-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 cisregulatory change in the Mdr65 gene [15]. In order to identify gene-regulatory changes on a whole-transcriptome scale in aamanitin-resistant D. melanogaster larvae, we performed a microarray study. As starting material, we used the most resistant of the four described a-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 (LC 50 = 2.16 mg/g of food) is at least as resistant to a-amanitin as the original Ama-KTT stock (LC 50 = 1.84 mg/g of food) (Figure 1), indicating that at least the majority of the aamanitin 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 aamanitin (LC 50 = 0.042 mg/g of food, data not shown).
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 a-amanitin. The complete set of differentially expressed genes (DEGs) can be found in Table S1. The isochromosome stock Ama-KTT/M/2 (LC 50 = 2.16 mg/g of food) was used as the experimental stock and has a 77.1 times higher LC 50 to a-amanitin than our sensitive control stock Canton-S (LC 50 = 0.028 mg/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 a-amanitin-containing food (at 1.5 mg/g of food, a concentration that is slightly below the LC 50 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 fullygrown 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.  . 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 a-amanitin, as symbolized by the yellow color. The larvae of group 3 (Ama-KTT/M/2) were treated with a-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. doi: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 a-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][19][20][21][22][23][24], while CYP6A2 metabolizes insecticides in enzyme assays [19,25]. Cyp12d1 is also associated with insecticide resistance [22,23,[26][27][28][29][30][31] and stress response [29][30][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][38][39], and Cyp305a1 (ecdysteroid synthesis and lipid storage regulation) [40].

Genes Encoding Cytochrome P450s and GSTs Are Inducible by a-Amanitin
Our next question was what genes are inducible by a-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 a-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.8fold) 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 a-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 a-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 a-amanitin in the resistant stock. As a result, nine genes were 2.1-to 3.8-fold downregulated in response to a-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 a-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.

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 a-amanitin in a dominant fashion [14]. Eighteen years later, a Californian D. melanogaster stock showed aamanitin 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 cisregulatory change in the Mdr65 gene is responsible for the resistance a-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 a-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.2fold down-regulated in response to a-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 upregulation 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 a-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 a-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 aamanitin (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 a-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 aamanitin. Mdr genes are thus not likely to confer the a-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 Type I and II detoxification, Mdr, and transcription factor genes with possible functions in detoxification processes are shown. The at least 2-fold differentially expressed genes are sorted by positive and negative fold-changes, followed by the genes that are not significantly differentially expressed. All p-values are corrected. The chromosomes, FlyBaseID, and probe ID numbers are presented. doi:10.1371/journal.pone.0093489.t001 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 a-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 a-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 a-amanitin.

Genome Enrichment Analysis Confirms 30-Year-Old QTL Mapping Data
In order to identify the regulatory pathway components that lead to the a-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 a-amanitin resistance to the polytene bands 95 and 98 on chromosome 3, we found signatures for both constitutive (group 2 versus group 1) and a-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 a-amanitin resistance is genetically controlled. The remaining four clusters of differentially expressed genes responded to a-amanitin in the larval food. The most interesting induced gene cluster is situated at cytological band 66A, which is close to Mdr65containing region 65A10 on the left arm of chromosome 3, to which a-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.8and 7.1-fold inducible by a-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 Table 2. Single gene analysis for Ama-KTT/M/2 on a-amanitin versus Ama-KTT/M/2 on no toxin (group 3 versus 2).
Type I and II detoxification, Mdr, and transcription factor genes with possible functions in detoxification processes are shown, sorted by positive and negative foldchanges. The at least 2-fold differentially expressed genes are sorted by positive and negative fold-changes, followed by the genes that are not significantly differentially expressed. All p-values are corrected. The chromosomes, FlyBaseID, and probe ID numbers are presented. doi:10.1371/journal.pone.0093489.t002 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 a-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.

Gene Ontology Enrichment Analysis Suggests Additional a-Amanitin Resistance Mechanisms
In order to explore if multiple mechanisms confer the resistance phenotype to a-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 a-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 a-amanitin. 2) 'Transferase activity' genes (GO 0016740) were on average 4.6fold higher expressed in Ama-KTT/M/2 (p = 7.61E-11), confirming our single gene analysis for the Gst and Ugt genes (Table 1) (Table S2). It is possible that cuticular proteins provide a protective layer against a-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 a-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 a-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 subprocesses. 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 a-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 upregulated, 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 aamanitin outside of tissues or perhaps even modifying it so that it becomes less toxic.
Next, we aimed to identify the gene ontologies that respond to a-amanitin in the resistant stock Ama-KTT/M/2. We thus compared Ama-KTT/M/2 on a-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 a-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 a-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 a-amanitin is a peptide, peptidases are good candidates to cleave it. To date, however, no specific enzyme is known that can inactivate a-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 a-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 a-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 a-amanitin and that yellow could play a protective role. Our gene ontology enrichment analysis further identified cellular components that respond to a-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 a-amanitin with a several thousand-fold induction of lipid particle genes suggests that cytoplasmic lipid particles contribute to the resistance to a-amanitin in the Ama-KTT/M/ 2 stock.

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 a-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, con-  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

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 a-amanitin to Ama-KTT/M/2 (group 3 versus group 2), Cyp316a1, 6d2, 4d8, 28d1, and 6t1 were upregulated 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.

Several Mechanisms Seem to Confer a-Amanitin Resistance
a-Amanitin is the principal toxin in some of the most deadly poisonous mushrooms, which inhibits the function of RNApolymerase 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 a-amanitin-inducible mechanisms that can explain the resistance to a-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 a-amanitin-inducible toxin entry blockage, mediated by cuticular proteins, the MRJP domain of the Yellow protein family, and Sgs proteins, 2) constitutive and a-amanitin-inducible phase I and II detoxification, mediated by the Cytochrome P450, GST, and UGT enzyme families (likely followed by excretion), 3) a-amanitin-inducible lipid particle gene induction, possibly leading to the sequestration of aamanitin in cytoplasmic lipid particles, and 4) a-amanitininducible peptidase genes, perhaps leading to the digestion of aamanitin either inside or outside (e.g. gut lumen) of cells (Figure 7).
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 a-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 a-amanitin-inducible salivary gland secretion genes Sgs1, Sgs3, Sgs5, Sgs7, and Sgs8 (each .300-fold induced) could perhaps bind to a-amanitin and reduce its uptake in the midgut. Another possibility is that a-amanitin is simply a stress factor that induces these and other genes. After all, a-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][54][55]. Cytochrome P450 proteins, which are encoded by Cyp genes, are known for their broad range  (Table 1). It is thus possible that one or all three of these genes contribute to the resistance to aamanitin. There is also evidence that Cyp12d1 is inducible by environmental stress factors, such as heat, oxidative stress, and air pollutants [29][30][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][60][61][62][63][64][65][66][67][68][69][70][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 a-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, 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 a-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. doi:10.1371/journal.pone.0093489.g006 however, possible that the a-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 a-amanitin in the larval food: the possible sequestration of aamanitin in lipid particles and the cleavage of a-amanitin by peptidases. A group of genes involved in the cellular component 'lipid particle' were on average more than 5200-times induced by a-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 a-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 a-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 a-amanitin resistance in D. melanogaster. In 1982 and 2000, two studies based on QTL mapping suggested that a-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 a-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 a-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 aamanitin 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 a-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 a-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 aamanitin. Furthermore, because we used whole larvae in our study, we cannot determine the relative importance that the different tissues play in the resistance to a-amanitin.
Our microarray data analysis did not reveal any gene-regulatory pathways that lead to the resistance to a-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 a-amanitin. However, the expression levels of both Hr96 and cnc revealed nothing that would lead us to conclude their role in a-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 a-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 a-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 aamanitin is mostly a constitutive trait.
Based on the mapping data from the two previous studies, we expected to find the a-amanitin resistance-conferring genes on Figure 7. A model of the four mechanisms that contribute to the resistance to a-amanitin in concert. The bicyclic octapeptide aamanitin is shown as a red 8. Cuticular proteins block some of the aamanitin from entering the cells (blockage). a-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. doi:10.1371/journal.pone.0093489.g007 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 a-amanitin in the Ama-KTT/M/2 isochromosome stock, future microarray studies should include samples of larvae that have been exposed to a-amanitin for different periods of time. Because first instar larvae are very small, the exposure to a-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 a-amanitin resistance in mycophagous Drosophila species, using the RNA sequencing approach. Mycophagous species are several orders of magnitude more resistant to aamanitin than D. melanogaster [11][12][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 a-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 a-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 a-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 b-amanitin, phalloidin, ibotenic acid, and muscimol should be used to test if crossresistance or independent mechanisms provide protection against the variety of mushroom toxins that mycophagous larvae encounter in their food source. Another pressing question is where a-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 a-amanitin could be a means to answer this question, but the analysis of the data might prove very difficult.

Conclusions
We suggest that the a-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 a-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 a-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 aamanitin. Future studies should include a time series of a-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.

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.  [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 a-Amanitin
In order to quantify and compare the levels of a-amanitin resistance of the D. melanogaster stocks, dose-response experiments were performed, which measured the survival from freshlyhatched first-instar larvae to adulthood. Flies able to completely hatch from their pupae were scored as survivors. The a-amanitinresistant stocks Ama-KTT and Ama-KTT/M/2 were tested on 11 a-amanitin concentrations, using 0 to 10 mg of a-amanitin per g of food in 1 mg increments. The a-amanitin-sensitive wild-type stocks Canton-S and the multi-balancer stock w  [3]; MKRS/TM6B, Tb [1] were initially tested on five concentrations ranging from 0 to 4 mg of a-amanitin per g of food in 1 mg 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 mg of aamanitin 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 25uC, 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 mL sterile Milli-Q water with or without dissolved a-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 LC 50 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 LC 50 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 a-amanitin-sensitive andresistant 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 a-amanitin, Ama-KTT/M/2 larvae were raised on 1.5 mg of a-amanitin per g of food (group 3), which is slightly lower than the LC 50 concentration of this stock. All larvae were raised until they reached the late third instar at 25uC, 70% humidity, and a 12:12 hour day/night cycle. Because not all larvae survived in the tubes and the larvae on a-amanitincontaining 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 highquality microarray data were used for normalization. Probesetlevel 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 a-amanitin slightly below the LC 50 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 a-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 aamanitin 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 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 a-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]. (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 a-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. (XLSX)