Transcriptome Analysis of Male Drosophila melanogaster Exposed to Ethylparaben Using Digital Gene Expression Profiling.

Ethylparaben (EP) has been shown to have estrogenic effects and can affect the normal development, longevity, and reproductive system of some animals. In this study, we investigated the effects of EP in male Drosophila melanogaster using transcriptome analysis or digital gene expression profiling. We then screened differentially expressed genes (DEGs) between the two groups (EP-treated and control group) of Drosophila, and performed clustering analysis, gene ontology (GO) function annotation, kyoto encyclopedia of gene and genomes metabolic pathway analysis. We found that EP enriched GO in three processes: cellular component, molecular function, and biological process. Consequently, we detected 13,959 genes and among them, 18 genes were identified to be significantly expressed between the EP-treated and control samples. Of these, seven genes were down-regulated, and eleven genes were up-regulated in EP-treated samples. Furthermore, four DEGs including two down-regulated genes (CG9465, CG9468) and two up-regulated genes (TotA, Sqz) were verified by real-time quantitative PCR. This study revealed the impact of EP on gene expression in fruit fly and provided new insight into the mechanisms of this response, which is helpful for understanding EP toxicity to humans.

Parabens, 4-hydroxy benzoic acid and alkyl or benzyl ester, have been used for about one century in cosmetics, food preservatives, and pharmaceutical and personal care products (PPCPs, Haman et al. 2015). Parabens include methylparaben (MP), ethylparaben (EP), propylparaben (PP), butylparaben (BP), etc., and have wide antimicrobial spectrum, high pH stability, and low cost (Boberg et al. 2010). Parabens are present with low concentration in environmental samples such as water (Stuart et al. 2012), air (Ramírez et al. 2010), soil (Nunez et al. 2008), and dust (Błedzka et al. 2014). In recent years, parabens are also detected in human. For example, it was detected in human breast tissue (Harvey and Everett 2004), milk (Schlumpf et al. 2010), and urine and serum (Carrasco-Correa et al. 2015).

Parabens are considered to be potential endocrine-disrupting chemicals (EDCs) by the Endocrine Society (Diamanti-Kandarakis et al. 2009). The measurement of the paraben is mostly in ng/liter, and paraben is considered to be endocrine disruptors (Haman et al. 2015, Petala et al. 2015). Some studies found paraben could promote the proliferation of human breast cancer cells MCF-7 (Charles and Darbre 2013, Khanna et al. 2014), affect the fecundity in male mice (Zhang et al. 2014), inhibit the mitochondrial respiratory capacities (Nakagawa and Moldéus 1998), and induce oxidative stress in the cells (Nishizawa et al. 2006, Pop et al. 2011, Shah and Verma 2011). Furthermore, parabens have estrogenic activity (Routledge et al. 1998, Darbre and Harvey 2008, Lange et al. 2014) and antiandrogenic activity (Satoh et al. 2005, Chen et al. 2007). Therefore, we hypothesized that parabens execute abovementioned function through influencing gene expression, just like steroid hormone.

To address this issue, we investigated the effect of paraben on gene expression in Drosophila melanogaster. We used D. melanogaster as the experimental organism for three reasons: (1) D. melanogaster can be used as a rapid, robust, reliable, and cost-efficient in vivo model for the possible toxic effects of parabens. (2) The use of Drosophila is ethically less controversial than that of other higher animal for in vivo studies. (3) The sequence of fruit flies is similar to human sequence and nearly 75% of human disease gene sequences are similar with D. melanogaster (Reiter et al. 2001, Lloyd and Taylor 2010), which supports its value for studying other human pathologies (Pandey and Nichols 2011, Prüßing et al. 2013, Gonzalez 2013, Alaraby et al. 2015). Our lab has been studying the toxicity of parabens in Drosophila in the past few years. We found that MP has significant toxic effects on the growth, development and fecundity in Drosophila, and the concentration of 200 mg/liter of MP has estrogenic effects in fruit flies (Gu et al. 2009). Different concentrations of PP can impact the reproduction and lifespan of Drosophila, and reduce the resistance to stress of heat shock and starvation in fruit flies (Li et al. 2015). PP can also affect mitochondrial functions, such as decreasing activities of superoxide dismutase (SOD) and catalase (CAT) and increasing malondialdehyde (MDA) content in Drosophila (Gu et al. 2010). And PP can significantly inhibit the gene expression of MnSOD and CuZnSOD (Jiang et al. 2011). EP have the reproductive toxicity on male Drosophila (Zhao et al. 2014), and can affect the gene expression of ERR, EcR, and YPR in D. melanogaster (Liu et al. 2014). The combination of MP and EP was found to affect preadult development period and lifespan of D. melanogaster (Chen et al. 2016). However, the mechanism behind parabens effect, especially the gene expression regulation, is still not clear. So we used digital gene expression (DGE) tag profiling to explore the mechanism by which the paraben affects D. melanogaster.

DGE technology can identify millions of differentially expressed genes (DEGs) before prior annotations and it allows the analysis for the organisms without much genetic information (Liu et al. 2014). Thus, DGE has been widely used to monitor the differences of transcriptional responses in different tissues and organs pressure in silkworm (Gao et al. 2013), cotton (Wei et al. 2013), moss (Nishiyama et al. 2012), Sagittaria trifolia (Cheng et al. 2013), Populus simonii and Populus nigra (Chen et al. 2012), Brassica napus (Jiang et al. 2013), Manduca sexta (Cao and Jiang 2015), Protaetia brevitarsis seulensis (Kyeongrin et al. 2015), and Cylas formicarius (Ma et al. 2016).

In this study, DEGs in male D. melanogaster in response to EP treatment were examined using DGE tag profiling technology. Then, we identified significantly expressed genes between EP treated and control groups. Finally, the results of DGE were validated in some up-regulated and down-regulated genes by qRT-PCR.

Materials and Methods

Drosophila Stocks and Treatment

In this study, we used the strain canton S of D. melanogaster. The fruit fly was reared on a standard corn meal, including corn flour, sugar, agar, yeast, and propionic acid, at 25°C, with a 50–60% humidity and 12:12 (L:D) h in the incubator.

When corn flour was prepared to feed Drosophila, EP was added directly to a final concentration 700 mg/liter for treated group, which is the maximum dose of food additives in China (). The control of Drosophila cultured in standard corn flour without EP.

We collected virgin flies within 8 h from the standard culture medium, and transferred them to the medium containing the EP (700 mg/liter) for 25 days. During the exposure, flies were transferred to the bottle with the same but new diet every two days to keep the food fresh. After 25 days, we isolated the EP-treated and nontreated (control) flies for subsequent experiments. In order to avoid individual differences, 25 fruit flies were randomly selected for each group.

RNA Isolation, Library Construction and Sequencing

We extracted the total RNA from flies using RNAiso plus kit (TAKARA, Dalian, China) according to the manufacture’s instructions. The fruit fly of each group was collected three times to extract the RNA. The integrity and quantity of RNA was inspected by 1.0% agarose gel electrophoresis in TAE and visualized under UV irradiation. The purity of RNA was examined by absorbance ratio between 260/280 nm in spectrophotometer (Daines et al. 2011).

The RNA from EP-treated and control group were used for library construction. Each group was repeated once. mRNA was enriched by magnetic beads and broken into short segments by adding the fragmentation buffer. By using random primers (random hexamers) of six bases the first strand cDNA can be synthesized. The second strand cDNA synthesis was performed by DNA polymerase I. Then the double-stranded cDNA was purified by AMPure XP beads. After end repair, poly A was connected to the segment and the fragments with specific size were selected by AMPure XP beads. Finally, a cDNA library was obtained by enrichment PCR and cleaned by PAGE gel electrophoresis. Each library preparation was sequenced on an Illumina Hiseq 2000 (Illumina, San Diego, CA; Chen et al. 2010).

Mapping Analysis and DGE Reads

For each DGE library, inferior label and adapter sequence tags were filtered. Then the clean tags were reserved and mapped to the reference genome by Top Hat 2. The algorithm of Top Hat 2 consists of three sections. First, the sequencing and transcriptome sequencing make a comparison. Second, the entire length of the sequence was compared with the exome of genome. Thirdly, the section of the sequencing was compared with the two exomes of genome (Kim et al. 2013).

Transcriptome data are often normalized for the length of the transcript as reads per bases per million reads (RPKM, Afrose et al. 2010). The results of genetic sequencing depth and length can be calculated by RPKM for reads count in the meantime. In general, it is the most common method of detection of gene expression (Mortazavi et al. 2008). In other studies, the balance of false positive and false negative number 0.3 RPKM establishing thresholds to determine whether or not a particular gene is expressed (Sam et al. 2011, Beyer et al. 2012, Wright et al. 2013). Differential expression analyzes between EP-treated and control groups were performed by the DEG Seq R package v1.12.0 (Anders and Huber 2010). The P-values were corrected by the Benjamini and Hochberg method (Benjamini and Hochberg 1995). Adjusted P-values are also called q-values. The threshold for the DGE is the q-values ≤ 0.05 and log 2 (fold-change) without limitations (Zhou et al. 2014).

Functional Analysis of DEGs

According to gene ontology (GO) terms, the differentially gene expressions were annotated in biological process, molecular function and cellular component (Huang et al. 2009). DEGs of GO enrichment analysis were carried out by GO Seq (Young et al. 2010), in which gene length and corrected deviation P-value <0.05 was deemed to have greatly enriched the differential gene expression. Kyoto encyclopedia of gene and genomes (KEGG) is a database resource used to understand biological system of the advanced functions. The formula was (Guan et al. 2013, Bai et al. 2014): where n is the number of DEGs of N, M is the number of genes with specific GO/KEGG pathway annotations, and m is the number of the DEGs with specific GO/KEGG pathway annotations of M. In this study, P-values ≦ 0.05 showed significantly pathways and enriched GO terms.

Quantitative Real-time PCR for Validating the DEGs

Real-time quantitative PCR (qRT PCR) analysis was used to verify the DGE results. Total RNA was the same for the DGE library. The single-strand cDNA was synthesized by PrimeScript RT Master Mix (TAKARA) in each group. The qRT-PCR was performed in CFX96 instrument (Bio-Rad, Hercules, CA). Specific oligonucleotide primers for some genes and ribosomal protein 49 (rp49) were designed by Primer Premier 5.0 software. The amplifications and measurements of qRT-RCR (Bustin et al. 2009, Bustin et al. 2010) by following the manufactures’ instructions were conducted using an AB7500 Real-Time PCR System (Applied BioSystems). Each reaction was run as following: denaturation program (95 °C for 5 min), followed by 35 cycles of 94°4 for 30 s, 60°0 for 34 s and 72°2 for 45 s with a primer pair in 20 µl reaction volumes. The program of melting curve was from the lowest 60°0 to the highest 95°5, with increasing 0.2°. per 15 s, to detect the identity and specificity of the qRT-PCR in the end stage of amplification. Each sample was technically performed in triplicate. The relative gene expression changes were calculated by the comparative CT (2-△△CT; Schmittgen and Livak 2008) and all samples were normalized to levels of rp49.

Statistical Analysis

All results were expressed as mean ± standard error (SE). All data were evaluated by the analysis of variance (One-way ANOVA) and SPSS software 18.0. Statistical significance level was considered significant at P < 0.05 or P < 0.01. The figures were draw by Origin8.0.

Results

Analysis of DGE Libraries

To establish a complete quantitative and qualitative gene expression response profile to EP, two DGE libraries were sequenced from two independent samples of male Drosophila with or without EP treatment. The main characteristics of these libraries are summarized in Table 1. The two EP-treated fruit flies’ libraries (termed EP-DM1 and EP-DM2) and two control one (termed Ctl-DM1 and Ctl-DM2) were checked their purification and quality. We have submitted our raw sequence data to NCBI, and the project accession number is PRJNA385568 or SRP106521. The Ctl groups showed 13,300,194, and 9,395,846 raw reads, of which 340,401 and 239,204 were multiple reads that were mapped to 12,037,198 and 8,515,079 unique tags, respectively; although the EP-treated groups showed 8,073,034 and 6,467,432 raw reads, of which 185,608 and 156,740 were multiple reads that were mapped to 7,324,138 and 5,875,615 unique tags, respectively.

Table 1.

Basic characteristics of tags in four libraries and sequencing reads mapping to the reference genome

Sample name	Ctl-DM1	Ctl-DM2	EP-DM1	EP-DM2
Total clean reads	13,236,139	9,351,457	8,059,462	6,456,485
Total mapped	12,377,599 (93.51%)	8,754,283 (93.61%)	7,509,746 (93.18%)	6,032,355 (93.43%)
Multiple mapped	340,401 (2.57%)	239,204 (2.56%)	185,608 (2.3%)	156,740 (2.43%)
Uniquely mapped	12,037,198 (90.94%)	8,515,079 (91.06%)	7,324,138 (90.88%)	5,875,615 (91%)
Reads map to “+” a	6,033,489 (45.58%)	4,263,653 (45.59%)	3,671,979 (45.56%)	2,942,992 (45.58%)
Reads map to “-” b	6,003,709 (45.36%)	4,251,426 (45.46%)	3,652,159 (45.32%)	2,932,623 (45.42%)
Nonsplice reads c	11,086,504 (83.76%)	7,828,127 (83.71%)	6,730,954 (83.52%)	5,396,636 (83.58%)
Splice reads d	950,694 (7.18%)	686,952 (7.35%)	593,184 (7.36%)	478,979 (7.42%)

Refers to sense strands.

Refers to anti-sense strands.

Refer to reads for the entire sequence is mapped to one exon.

Splice reads also called junction reads, refers to reads mapped to the border of exon.

Analysis of Gene Differential Expression

For analysis of gene expression, RPKM can be calculated and standardized cleaning clear tag number of each gene. RPKM distribution of all genes and a RNA-seq correlation test were examined the difference of EP-treated samples and Ctl samples. In EP-treated and Ctl samples, the Pearson’s correlation coefficients were 0.979 and 0.985, respectively. RPKM density distribution showed a number of different genes in EP-treated and Ctl samples.

These results showed DEGs in Ctl and EP-treated libraries. The transcriptome tested with log 2 (fold-change) without limitations and q-value < 0.05 was used to compare the experimental library to the control one. Consequently, 18 significantly expressed genes were identified from 13,959 genes screened between the EP-treated and Ctl samples (P-values < 0.05). Among them, seven genes were down-regulated, and eleven genes were up-regulated in EP-treated samples, which could be found in Table 2.

Table 2.

List of all the different genes between control and ethylparaben group and pathways in which the different genes related to EP exposure are active

Gene ID	Gene name	Readcount EP_DM	Readcount DM	Log 2 fold change	Padj	Pval	Regulated	Gene annotation
FBgn0010768	squeeze	114.41	49.19	1.22	0.03	3.35E-05	Up	Metal ion binding, nucleic acid binding, sequence-specific DNA binding, transcription factor activity.
FBgn0028396	Turantot A	1,194.98	684.24	0.80	0.01	6.33E-06	Up	Elemental activities, such as catalysis or binding, describing the actions of a gene product at the molecular level.A given gene product may exhibit one or more molecular functions.
FBgn0032067	CG9465	9.56	83.35	−3.12	0.01	9.71E-07	Down	Alpha-mannosidase activity, carbohydrate binding, zinc ion binding.
FBgn0032068	CG9466	914.92	2,017.87	−1.14	0.01	2.04E-07	Down	Alpha-mannosidase activity, carbohydrate binding, zinc ion binding.
FBgn0032069	CG9468	1,675.59	4,500.58	−1.43	2.01E-13	1.44E-17	Down	Alpha-mannosidase activity, carbohydrate binding, zinc ion binding.
FBgn0032144	CG17633	1,471.66	862.90	0.77	0.01	1.10E-05	Up	Metallocarboxypeptidase activity, zinc ion binding.
FBgn0032147	IP3K1	473.77	266.86	0.83	0.04	5.81E-05	Up	NOT calmodulin binding, inositol-1,4,5-trisphosphate 3-kinase activity.
FBgn0033294	Mal-A4	273.04	1,097.74	−2.01	0.01	8.20E-06	Down	Alpha-glucosidase activity, cation binding.
FBgn0033592	CG13215	63.08	16.19	1.96	0.01	1.16E-06	Up	Protein coding.
FBgn0035154	CG3344	409.60	218.41	0.91	0.01	4.05E-06	Up	Serine-type carboxypeptidase activity.
FBgn0035358	CG14949	534.88	275.30	0.96	0.04	5.53E-05	Up	Protein coding.
FBgn0036619	Cpr72Ec	330.99	102.26	1.69	2.94E-05	6.33E-09	Up	Structural constituent of chitin-based cuticle.
FBgn0036764	CG5535	314.55	163.56	0.94	0.02	2.01E-05	Up	Amino acid transmembrane transporter activity, cationic amino acid transmembrane transporter activity.
FBgn0037782	Npc2d	6.80	88.36	−3.70	1.73E-11	2.48E-15	Down	Sterol binding.
FBgn0051233	CG31233	2,936.77	1,808.36	0.70	0.03	2.74E-05	Up	Aminopeptidase activity, metallopeptidase activity, zinc ion binding.
FBgn0051664	CG31664	35.15	102.70	−1.55	0.01	1.15E-05	Down	Protein coding.
FBgn0085232	CG34203	76.97	27.15	1.50	0.02	2.27E-05	Up	Protein coding.
FBgn0263336	CR43417	257.68	597.55	−1.21	0.01	1.17E-06	Down	ncRNA, putative noncoding RNA.

Functional Annotation of DEGs

The ontology database covers three domains: cellular component, molecular function, and biological process. The DEGs were checked by GO analysis and there were differently enriched GO terms in molecular function or biological process (Fig. 1). For example, there were three, six, and fourteen genes in the biological process, cellular component and molecular function enriched GO terms, respectively (Fig. 2). The results showed that gene differentially expressed between the EP-treated and Ctl group, and some genes enriched GO in cellular component, molecular function and biological process participated EP response (Table 2), such as sterol transport, response to oxidative stress, mannosidase activity, mannose metabolic process, metal ion binding, zinc ion binding and transition metal ion binding.

Fig. 1.

GO functional enrichment of genes differentially expressed between the control and EP-treated group in biological process (a), cellular component (b), molecular function (c). Gene classification based on gene ontology (GO) for differentially expressed genes in the two groups. The number of genes in GO terms was analyzed using GO Slim Assignment.

Fig. 2.

GO analysis of differentially expressed genes between the control and EP-treated group. The differentially expressed genes are classified into three categories, cellular component, molecular function, and biological process. The flowing axis shows number of genes, with ‘*’ for the enrichment of GO term.

KEGG pathway can analyze differently expressed genes. The classifications showed that some genes expression intensity changed in relation with sucrose metabolism, phosphatidylinositol signaling system and other pathways.

DGE Results Were Confirmed by qRT-PCR

In order to validate DGE gene results, four genes were selected from 18 DEGs, including two up-regulated (TotA, Sqz) and two down-regulated (CG9465, CG9468), for qRT-PCR in the EP-treated and Ctl group and the sequences are listed in Table 3. The expression patterns of the four genes were consistent with the tag-mapped pattern (Fig. 3), which increased the credibility of the DGE in this study.

Table 3.

The primers used in qRT-PCR

Gene	Primer sequence
rp49-F	5′-ATGACCATCCGCCCAGCATAC-3′
rp49-R	5′-GCATCAGATACTGTCCCTTGAAGC-3′
ToA-F	5′-TCAATAGCACCCAG GAACTAC-3′
ToA-R	5′-AGAATCCAGTTGCAAGGC-3′
Sqz-F	5′-CGAGC AGCAGGAACAGAAGT-3′
Sqz-R	5′-TCGCAAAACAATCCGAAGCA-3′
CG6465-F	5′-TCGCAAAACAATCCGAAGCA-3′
CG9465-R	5′-TTGTTTC GGGACATGACTGC-3′
CG9468-F	5′-ATTTGTGCGGCATTGTGCTT-3′
CG9468-R	5′-TCCTGTCGCATTCTGGCATT-3′

Fig. 3.

qRT-PCR validations of expression levels of candidate genes from DGE analysis. Quantitative RT-PCR validation of differentially expressed genes in control (a basal diet) and treatment (diet with 700 mg/liter EP) samples. Up-regulated genes (a): ToA; (b): Sqz; down-regulated genes; (c): CG9465; (d): CG9468). All data were normalized to the expression level of actin. Data represent fold change of relative quantification in treatment and control groups. The error bars represent the range of the fold change as determined by the data assist software.

Discussion

In this study, the male Drosophila was exposed to EP and transcriptome analysis was performed by using DGE profiling. EP belongs to parabens, and acts as EDCs. Further research is needed to provide new insights into the mechanism(s) through which endocrine disruptors elicit their effects on biological systems and human and animal health (Vo et al. 2010). We investigated the influence of EP on Drosophila. EP have the reproductive toxicity on male Drosophila (Zhao et al. 2014), and can affect the gene expression of ERR, EcR, and YPR in D. melanogaster (Liu et al. 2014). However, the mechanism by which EP affect the genes expression was still unknown. Considering the molecular mechanism of these processes, transcriptome sequencing was performed and compared for each sample. On the basis of GO enrichment analysis, there were many DEGs that enriched GO in molecular function, cellular components and biological processes. And they were mainly shown in mannose metabolic process, lytic vacuole, lysosome, alpha-mannosidase activity, mannosidase activity, hydrolase activity, hydrolyzing O-glycosyl compounds, hydrolase activity, and acting on glycosyl bonds. We selected the TotA, Sqz, CG9465, CG9468 and the expression pattern was further confirmed by qRT-PCR. Thus, we identified the candidates for genes involved in the EP treatment.

The previous studies showed that Turandot A (TotA) is one of the eight members belonging to the Tot genes family distributed at three different sites in the Drosophila genome. All Tot genes are under pressure, such as heat shock, bacterial infection, exposure to ultraviolet light, or paraquat feeding, suggesting that all members of this family have played a compressive resistance in Drosophila (Petre et al. 2012). TotA is a member of a family of Drosophila humoral compressive resistance genes which hasn’t be described yet in mammalian homologs (Ekengren and Hultmark 2001). TotA is composed of a secreted protein of 129 amino acids and is well known to respond to environmental pressures such as chemical insult, heat shock, bacterial infection and ultraviolet irradiation (Ekengren et al. 2001). TotA has been the best explanation to the conservation of its expression in response to cellular stress pathways (Agaisse and Perrimon 2004). Some studies suggested that in the development process of larvae to adult flies, the TotA gene has a higher expression in the tolerance as compared with the nontolerance (Mahapatra and Rand 2012). However, our laboratory has found that mixing the MP with the EP can affect the larval to adult developmental stages (Chen et al. 2016). Several pathways are associated with TotA gene expression. Mekk1 (Mitogen-activated protein kinase kinase 1), an integral part of the MAPK pathway, is considered to be the regulated genes of TotA upon septic injury in Drosophila (Brun et al. 2006). In Drosophila, Imd and JAK–STAT pathways upon septic injury can also regulate TotA expression (Agaisse et al. 2003). The JAK–STAT pathway can be affected by many environmental poisons, such as cadmium (Monroe and Halvorsen 2009). TotA as an indicator will have a wider impact mechanism to address environmental chemical attack in future analysis (Mahapatra and Rand 2012). Our study, suggests TotA, an up regulated gene in EP-treated group, works as a regulatory cytokine gene in the role of conservation/immune signaling pathway to regulate EP toxicity.

As far as down-regulated genes were concerned, CG9465 and CG9468 can influence alpha-mannosidase activity, carbohydrate and zinc ion binding. The additional genes subjected to glucose repression (CG9468), the actual pool of affected carbohydrases and lipases could be potentially larger (Chng et al. 2014). However, dietary protein:carbohydrate (P:C) balance is the main driving force of lifespan and the aging processes (Lee 2015). High concentrations of PP can reduce the lifespan although 1,000 mg/liter PP can extend the life in fruit flies (Li et al. 2015). The combination MP with EP can reduce the longevity of flies (Chen et al. 2016). We performed qRT-PCR to confirm the effectiveness of the DGE and obtained the final results. Therefore, we confirmed our sequencing results of authenticity.

We identified two KEGG pathways, metabolic pathways and signal transduction, that were greatly enriched, and thus, obviously active in EP-treated Drosophila. The first one involved in the operation of nuclear receptor steroid hormone pathway is recognized as the targets of EDCs in vertebrates (Maglich et al. 2001). Bisphenol A (BPA), similar with EP, has estrogen effect and can affect the development and growth of Drosophila by interfering with the endocrine signals and estrogen receptor gene expression in D. melanogaster (Weiner et al. 2014).

In short, we have described genes that are expressed differentially in the male fruit flies exposed to EP. And 18 different genes are related to cellular component, molecular function, and biological process in D. melanogaster. We acknowledge that the EP-response mechanism is rather complicated, but this is only a preliminary description. Further, in-depth study is indispensable to illuminate the function of 18 DEGs. The data provided by this study will be useful to understand the toxicity of EP in D. melanogaster from transcriptomic view.