Hydroxyl-substituted double Schiff-base condensed 4-piperidone/cyclohexanones as potential anticancer agents with biological evaluation

Article abstract of DOI:10.1080/14756366.2018.1501042

Novel hydroxyl-substituted double Schiff-base 4-piperidone/cyclohexanone derivatives, 3a–e, 4a–e, 5a–d, and 6a–c, were synthesized and fully characterized by ¹H NMR, IR and elemental analysis. The cytotoxicity against human carcinoma cell lines

Full text of DOI:10.1080/14756366.2018.1501042

Science of the Total Environment 661 (2019) 226–234
Contents lists available at ScienceDirect
Science of the Total Environment
journal homepage: www.elsevier.com/locate/scitotenv
Short Communication
Comparative transcriptomic analysis provides insights into the response
to the benzo(a)pyrene stress in aquatic firefly (Luciola leii)
a,
a
a
a
b
a,
Qi-Lin Zhang ⁎, Jun Guo , Xian-Yu Deng , Feng Wang , Jun-Yuan Chen , Lian-Bing Lin ⁎
a
Faculty of Life Science and Technology, Kunming University of Science and Technology, Kunming 650500, China
LPS, Nanjing Institute of Geology and Paleontology, Chinese Academy of Sciences (CAS), Nanjing 210008, China
b
H I G H L I G H T S
G R A P H I C A L A B S T R A C T
•
•
29,297 of 54,282 genes obtained in
Luciola leii transcriptome were anno-
tated.
Most DEGs were involved in xenobiotic
biodegradation, biomacromolecule me-
tabolism.
•
•
•
This study discovered a set of useful SSR
and SNP molecular markers.
Eight DEGs were identified as key BaP-
responsive genes.
This study firstly explored transcrip-
tional response to BaP in aquatic insects.
a r t i c l e i n f o
a b s t r a c t
Article history:
Many studies have reported that behavior and bioluminescence of fireflies could be affected by changes in envi-
ronment conditions. However, little is known about how the deterioration of the aquatic environment affects
aquatic fireflies, particularly with respect to molecular responses following exposure to water pollutants, such
as benzo(a)pyrene (BaP), which is a key indicator in environmental risk assessment because of the hazards it
poses. Here, whole transcriptome sequencing and gene expression analysis were performed on freshwater fire-
flies (Luciola leii) exposed to BaP (concentration of 0.01 mg/L). Four transcriptomic libraries were constructed for
the control and treatment groups, including two biological replicates. From the mixed pools (each pool contains
60 individuals from three time points), a total of 54,282 unigenes were assembled. Furthermore, 329,337 of
Single-nucleotide Polymorphisms (SNPs) and 1324 of Simple Sequence Repeats (SSRs) were predicted using bio-
informatics, which is useful for the future development of molecular markers. Subsequently, 2414 differently
expressed genes (DEGs) were identified in response to BaP stress in comparison to the control, including 1350
up-regulated and 1064 down-regulated DEGs. Functional enrichment showed that these DEGs are primarily re-
lated to innate immunity; xenobiotic biodegradation and response, biomacromolecule metabolism, biosynthesis,
and absorption. Eight key BaP-responsive DEGs were screened to survey the dynamic changes of expression in
response to BaP stress at different time points, and to validate the RNA sequencing data using quantitative
real-time PCR. The results indicate that the expression of genes encoding UGT, CYP3A, CYP9, CYP6AS5 and
ADHP were induced, while those encoding UGT2B10L, PTGDS, and ALDH were reduced, to participate in response
to the BaP exposure and potentially help counteract the adverse effects of BaP. This investigation provides insight
into the toxicological response of fireflies to the occurrence of water deterioration.
Received 13 December 2018
Received in revised form 7 January 2019
Accepted 13 January 2019
Available online 15 January 2019
Editor: Henner Hollert
Keywords:
Luciola leii larvae
Transcriptomics
Gene expression
Benzo(a)pyrene
Molecular toxicology
RNA sequencing
©
2019 Elsevier B.V. All rights reserved.
⁎
Corresponding authors.
E-mail addresses: zhangqilin88888@126.com (Q.-L. Zhang), linlb@kmust.edu.cn (L.-B. Lin).
https://doi.org/10.1016/j.scitotenv.2019.01.156
048-9697/© 2019 Elsevier B.V. All rights reserved.
0

Q.-L. Zhang et al. / Science of the Total Environment 661 (2019) 226–234
227
1
. Introduction
mechanisms of aquatic insects under water contaminant stress. Fur-
thermore, the genetic resources obtained in this study will be useful
for development of novel molecular markers and assessment of envi-
ronmental risk in water bodies.
The family Lampyridae (Insecta: Coleoptera: Cantharoidae), also
known as fire flies or bioluminescent beetles, consists of over 100 gen-
era and 2000 nominated species (Fu, 2005). Fireflies can be divided
into terrestrial, aquatic and semi-aquatic lineages according to the living
habits of their larvae. Most firefly larvae are terrestrial while a small
number of larvae live in freshwater (namely the aquatic fireflies); how-
ever, their survival strongly depends on the quality of the water body
2. Materials and methods
2.1. L. leii larvae samples and BaP exposure
(
Fu et al., 2005). Aquatic fireflies are very rare and are distributed only
Healthy 5th instar L. leii larvae (body size: 1.2 ± 0.2 cm) (Fig. S1)
were purchased from the Culture Centre of Fireflies, Ganzhou, Jiangxi
Province, China, and were maintained at the Evo-dev Institute of Nan-
jing University, Beihai, Guangxi, China, as previously described (Fu,
2005). The experimental larvae were separated into two jars
representing BaP-treated and control groups. A total of 100 individuals
were maintained under water in 1.5 L wide-mouth plastic jars contain-
ing 1.0 L sterilized lake water. BaP (B1760, Sigma-Aldrich, USA, purity
N99%; molecular weight: 252.31) was dissolved into dimethyl sulfoxide
(DMSO, Sigma-Aldrich, USA, N99% purity) to reach a stock concentration
of 200 μM. Subsequently, the BaP-treated groups were exposed to fresh-
water containing 0.01 mg/L of BaP. The group treated with an equal con-
centration (0.03%) of DMSO in freshwater was used as a control. The
optimal BaP concentration used for L. leii exposure was determined
based on previous studies (unpublished data). In brief, the concentra-
tion exhibiting the maximum number of BaP-effected-expressed cyto-
chrome P450s (CYP450s), which is a suitable molecular indicator to
evaluate BaP pollutants (Camatini et al., 1998; Ryeo-Ok et al., 2013),
was selected as the suitable concentration in this study. At 6, 12 and
24 h post-exposure (hpe), 20 L. leii larvae individuals were collected
from each of the treatment and control groups. The experiments were
independently repeated two times, representing two biological repli-
cates. All samples (whole bodies of L. leii larvae) were immediately fro-
zen in liquid nitrogen and stored at −80 °C for further RNA extraction.
in Asia and Jamaica. Thus far, only seven species of aquatic fireflies (all
of them belonging to the genus Luciola) have been described in detail
(
Luciolinae) is a widely investigated aquatic species of firefly that has
completely aquatic larval stages. L. leii is primarily distributed in the
middle and lower reaches of the Changjiang River in mainland China,
and it is used as an indicator of good water quality (Fu and Ballantyne,
Fu et al., 2004; Fu et al., 2010). Among these, Luciola leii (subfamily
2
006; Vongsangnak et al., 2016). However, the mechanism underlying
the molecular toxicology in response to stress caused by contaminants
in water bodies remains largely unknown. Moreover, the genes and
pathways potentially involved in this toxicological and biological re-
sponse still remain to be identified.
Polycyclic aromatic hydrocarbons (PAHs) are a diverse group of hy-
drocarbons with over 1000 different compounds that are composed of
two or more fused benzene rings (Abdel-Shafy and Mansour, 2016).
PAHs can enter water bodies through different routes and are found as
a mixture of various PAH compounds at differing levels (Xiao et al.,
2
018). Due to the adverse effects of PAHs on various organisms, includ-
ing humans, PAHs have become important environmental contami-
nants of concern. Many investigations demonstrate that PAHs are a
crucial driving factor in carcinogenic, mutagenic and teratogenic pro-
cess (Xiao et al., 2018). Recently, due to ecological deterioration and in-
tensifying human activity, the PAHs can easily be found in rivers, lakes
and oceans. Benzo(a)pyrene (BaP) is a PAH known as one of the most
potent environmental carcinogens, and human activity is the main
cause of its release into the aquatic environment (Xiao et al., 2018).
One survey of aquatic sediments in China showed that the BaP content
was 0.06–52.12 ng/g in the Yangtze River, 0.19–36.71 ng/g in the Liaohe
River, 0.7–69.4 ng/g in Taihu Lake and 18.18–701.12 ng/g in the Pearl
River Delta (Jia et al., 2012). The average concentration of BaP in the sur-
face runoff of Hangzhou city, China, was up to 1582 ng/L, and reached
2.2. RNA isolation, sequencing and transcriptome assembly
Total RNA was extracted from the samples (whole bodies) using the
Ambion RNA Pure Kit (Invitrogen, USA). Residual genomic DNA was re-
moved by RNase-free DNase (Qiagen, Germany). An Agilent 2100
Bioanalyzer (Agilent Technologies, USA) was used to confirm RNA
structural quality (RNA integrity number (RIN value) ≥ 7). The RNA
from 20 individuals that were collected at each time point was diluted
to the same concentration (300 ng/μL) with RNase-free water, and
these were then pooled equally into the control and treatment samples
that were used for further library construction. Library construction was
performed at the Beijing Genome Institute (BGI, China) using the
Illumina RNA Sample Preparation Kit (Illumina, USA) according to the
manufacturer's protocol. The four libraries were sequenced on the
Illumina HiSeqTM 2000 platform to obtain 100 bp paired-end reads
(BGI, China).
The data filtering and de novo assembly of raw data were imple-
mented under the default parameters unless otherwise stated. First,
the raw data was filtered by removing the adaptor sequences, reads
with N10% unknown bases (N) and low-quality reads using a FastQC
software (v 0.11.7), which was used during quality control to check
raw sequence data (Li et al., 2015). Next, de novo assembly of mixed
read pooling generated by a combination of the four libraries was per-
formed using Trinity (Grabherr et al., 2011; Haas et al., 2013). The
TGICL pipeline was then used to remove redundancy of assembled
unigenes by a fast clustering (Pertea et al., 2003). To perform annota-
tion, all unigenes were searched in non-redundancy (NR), gene ontol-
ogy (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and
Swiss-Prot databases using the BLASTX tool (threshold of 1.0e-5). To
evaluate the completeness of the assembly, the generated unigenes
3
00–400 ng/L in the source water of Meiliang Bay in Taihu Lake (Jia
et al., 2012). In addition, few of the currently available studies have in-
vestigated the metabolic effects of BaP on freshwater insects (Borchert
et al., 1997). However, several recent publications have explored
large-scale gene expression changes in response to BaP exposure
using RNA sequencing. These studies found that the genes expressed
in response to BaP exposure were primarily involved in drug metabo-
lism, xenobiotic metabolism and innate immune signaling in the
staghorn coral species Acropora hyacinthus embryo and the adult Chi-
nese amphioxus Branchiostoma belcheri (Xiao et al., 2018; Zhang et al.,
2
019a). To date, the molecular mechanisms involving toxicological
and biological responses to water contaminants is completely unknown
for aquatic insects. Thus, it is necessary to employ suitable insect
models, such as the aquatic firefly, to explore molecular toxicological
mechanisms in response to important water contaminants. Moreover,
the aquatic insect models and their genetic resources will be useful for
environmental monitoring and development of technology used for
assessing environmental risk in water bodies.
The goal of this study was to explore the transcriptional response of
L. leii larvae to BaP exposure at a whole-transcriptome level, as well as
identify BaP-responsive genes and reveal their functional information.
RNA sequencing technology was used to obtain transcriptomes of L.
leii larvae, and then mass gene expression profiling was analyzed in
combination with bioinformatics and quantitative real-time PCR (qRT-
PCR). The results will be helpful for understanding the toxicological
were compared with
a core set of arthropods using BUSCO
(Benchmarking Universal Single-Copy Orthologs) v3 (http://busco.

228
Q.-L. Zhang et al. / Science of the Total Environment 661 (2019) 226–234
ezlab.org/), which is based on evolutionarily informed expectations of
gene content, following previously published methods (Zhang et al.,
2.6. Validation of RNA-Seq data and observation of gene expression by qRT-
PCR
2
017b). The completeness of the unigene set was automatically evalu-
ated, and percentage of “complete”, “fragmented”, and “missing” were
obtained by automatic output of BUSCO software.
To validate the RNA-Seq data, 12 DEGs (eight key BaP-responsive
and four randomly selected DEGs) were used for qRT-PCR analysis
from a sample of the same RNA that was used for the RNA-Seq experi-
ments. Furthermore, the expression level of eight key BaP-responsive
DEGs were detected in a replicated experiment using 5th instar L. Leii
larvae exposed to BaP for 6, 12, 24 and 48 h (hpe). The sampling and
RNA purification was conducted following methods described in
Sections 2.1 and 2.2. Specific qRT-PCR primers were designed using Bea-
con Designer 7 (Tables S1). Synthesis of the first-strand cDNA and qRT-
PCR analysis was performed according to methods used in our previous
studies (Zhang et al., 2017e). The elongation factor 1A (EF1A) gene was
used as a housekeeping gene (Tang et al., 2017) and each qRT-PCR reac-
tion was conducted with three technical replicates and two biological
replicates. Data was visualized with the SDS tool (version 2.0) in the
ABI 7300 Real-Time PCR system. Gene expression normalization and
2
.3. Identification and analysis of SSRs and SNPs
Simple sequence repeats (SSRs), also known as microsatellites,
have been widely used as genetic markers, genetic evolution, gene
mapping, comparative genomics and population genetics in insects
(
Wang et al., 2009). In brief, MIcro SAtellite (MISA) (v1.0) software
was used to detect SSRs in all the assembled unigenes (Thiel et al.,
003). The SSR primer design was then performed using the Primer3
tool (Triinu and Maido, 2007). Single nucleotide polymorphisms
SNPs) are single nucleotide mutations that occur between samples
2
(
or individuals at the DNA or RNA level. Because of variations in
some key functional nucleotide mutations, SNPs may lead to changes
or losses of protein function, which can result in adaptive evolution
of insects to their habitats (Jha et al., 2016). Therefore, SNPs are use-
ful molecular markers during investigations of environmental adap-
tation using aquatic fireflies. Alignment of multiple clean reads
corresponding to unigenes was conducted using HISAT (v0.1.6-
beta) (Daehwan et al., 2015) and detection and filtering of low-
quality SNPs was conducted using GATK analysis (v3.4-0)
−
△△Ct
statistics were conducted based on the 2
method (Livak and
Schmittgen, 2001). Namely, Ct values of all technical replicates of each
gene in each sample were generated by PCR system. Then, in each bio-
−
△△Ct
logical replicate, mean value of Ct, △Ct, △△Ct and 2
values of
each gene were calculated. Next, mean and standard deviation (SD) of
−△△Ct
2
values from two biological replicates were obtained. The signif-
icance of differences between control and BaP-treated firefly at each
time point (6, 12, 24 and 48 hpe) was evaluated with one-way
ANOVA followed by Bonferroni post-tests using IBM SPSS statistics 22.
(
McKenna et al., 2010).
2
.4. Analysis of differentially expressed genes
3. Results and discussion
To quantify gene expression levels, clean reads of each library
3.1. Overview of RNA sequencing data
were first mapped to the reference unigene set of the L. leii tran-
scriptome assembled in this study using Bowtie 2 under the default
parameters (Langmead et al., 2009). Furthermore, the RSEM (RNA-
Seq by Expectation Maximization) tool was used to calculate the
fragments per kilobase of unigenes per million fragments mapped
A flow diagram of the experiments and bioinformatic analyses is
presented in Fig. 1A. After RNA sequencing and quality control, approx-
imately 3 Gb of clean reads were generated from each of the four librar-
ies, and the Q20 percentage of the clean reads in each library was higher
than 97.45% (Table 1). A total of 294,391,102 clean reads, with an aver-
age GC content of 41.24%, were generated from the four libraries
(Table 2). These clean reads have been submitted to the Sequence
Read Archive (SRA) database of NCBI (Accession numbers
SRR8145632-SRR8145635). A total of 54,282 unigenes and 105,644
transcripts was obtained, with an average unigene length of 999 bp
and an N50 of 2409 bp. Functional annotation analysis showed that
53.97% of the unigenes (29,297) were annotated to at least one database
of NR, Swiss-Prot, GO and KEGG, while the remaining unigenes were not
annotated. The length distribution showed that 16.43% of the unigenes
were N2000 bp (Fig. S2). Overall, the large-scale mixed assembly of
the four libraries generated a higher quality transcriptome than that
had been previously obtained by other studies of aquatic fireflies
(
FPKM) (Zhang et al., 2017d) according to the amount of mapping
reads on each unigene. Sample reproducibility was assessed be-
tween the two biological replicates based on the FPKM values of
each unigene. Differentially expressed genes (DEGs) were deter-
mined and extracted using DESeq software, which is an R package
used for identification of DEGs for transcriptomes containing biolog-
ical replicates (Anders and Huber, 2010). The DEGs identified by the
DESeq tool shows more conservative and lower false rates than those
identified by other software with similar function (e.g. DEGseq and
NOISeq) (Behringer et al., 2015). The threshold for DEG determina-
tion was set as fold changes (FC) ≥ 2 (|log2 ratio| ≥ 1) and 0.01 p-
values (Wald test in DESeq) corrected by the Benjamini-Hochberg
method (false discovery rate, FDR).
(
Luciola aquatilis) (Vongsangnak et al., 2016). For example, the N50
and mean values of unigene length in this study (2409 bp and 999 bp,
respectively) are longer than that (1889 bp and 978 bp, respectively)
of Vongsangnak et al. (2016); moreover, in this study (16.43%), the per-
centage of unigenes with N2000 bp is higher than that in previous pub-
lications (13.58%) (Vongsangnak et al., 2016). In addition, results of the
BUSCO analysis based on 2675 near-universal single-copy orthologs in
the arthropod gene set showed that for the L. leii transcriptome assem-
bled in this study, 68.4% of the unigenes were “complete”, 16.2% were
“fragmented” and the remaining genes were “missing”. In comparison
to the 46 other arthropod species assemblies (Theissinger et al., 2016),
the results show that the integrity of the assembled unigenes is higher
than that of the majority of assembled transcriptomes. Therefore, a
high-quality unigene set for L. leii larvae was obtained here and subse-
quently used for downstream analyses. The high quality unigene set
for L. leii larvae is a reasonable and novel outcome of the present work.
2
.5. Functional enrichment analysis of differentially expressed genes
To understand the main biological, cellular and molecular functions,
as well as the signaling pathways of DEGs, the corresponding GO and
KEGG annotations were extracted for each annotated DEG using custom
Perl scripts. Next, GO enrichment analysis of DEGs was performed using
the Blast2GO pipeline (Conesa et al., 2005) with Fisher's exact test (p-
value calculation), and they were enriched into three subcategories,
namely biological processes, molecular function, and cellular compo-
nents. Next, KOBAS 2.0 was used with the default options (http://
kobas.cbi.pku.edu.cn/) to conduct KEGG enrichment analysis for the
DEGs. The Benjamini-Hochberg adjustment was used for correction of
the significant level of all the GO and KEGG terms. The FDR threshold
was set as 0.05.

Q.-L. Zhang et al. / Science of the Total Environment 661 (2019) 226–234
229
A
B
Control
BaP treatment
6
12 24 hpe
6
12 24 hpe
RNA extraction
of 20 individuals
at each time point
RNA pooling
RNA pooling
Control
RNA pool
Treatment
Log10(FPKM) of control_1
C
Rep_1
Rep_2
Rep_1
Rep_2
De novo assembly
Annotation
Transcriptome
NR
Swiss-Prot
GO
Mapping
KEGG
Quantification of gene expression
Screening of differentially expressed genes
Target genes
GO enrichment
KEGG pathway
qRT-PCR
Log10(FPKM) of treatment_1
Fig. 1. (A) A flow chart of experiments and bioinformatic analyses. Hpe indicates hours post-exposure. (B) Correlation analysis between two control samples based on gene expression
levels. (C) Correlation analysis between two treatment samples based on gene expression levels. When black dots are close to the pink oblique line, this indicates that the two samples
are highly similar. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
3
.2. SSRs and SNPs identification
33.70% and 33.56%, respectively, while the other four SNP types (A/C,
A/T, C/G and G/T) represented 7.54%, 12.74%, 4.93% and 7.53%, respec-
tively. This study presents the first detection of molecular markers for
aquatic fireflies. Genetic diversity evaluation based on SSR and SNPs is
an effective approach for conservation biology research of aquatic ani-
mals (Miller et al., 2014). Therefore, these molecular markers are valu-
able for the protection of L. leii larvae in water environments.
In total, 1342 SSR markers were identified in 1278 unigenes. Among
these unigenes, 54 contained N1 SSR. Among those SSRs, the most fre-
quent motif type was tri-nucleotide (642, 47.84%), followed by di-
nucleotide (400, 29.81%), mono-nucleotide (255, 19.00%), quad-
nucleotide (26, 1.94%) and other (19, 1.42%) nucleotide motif types
(Fig. 2A). Moreover, the size of the SSR sequences was mainly distrib-
uted in the 12–24 bp range, occupying 99.25% of the total SSRs.
Among the SSRs, a total of 29 motif types were detected. The repeated
AAT/ATT nucleotide type was the most frequent motif type in this inves-
tigation (286, 21.31%), followed by A/T (253, 18.85%) and AT/AT (238,
3.3. Analysis of differentially expressed genes
Clean reads were mapped to the reference unigene set for each li-
brary and the mapping statistics are summarized in Table 1. The per-
centage of total mapped reads ranged from 83.46% to 88.51% for the
four libraries, of which approximately 65% and 22% accounted for
unique and multi-position matches, respectively, for each library. A
1
7.73%). In the current investigation, a total of 329,337SNP sites were
detected. The number of transition and transversion types was
21,500 and 107,837, respectively (Fig. 2B). The frequency of SNPs in
2
the L. leii larvae transcriptome was 1/165 bp, meaning that on average
an SNP site could be found per 165 bp of unigene length. The number
of A/G and C/T was maximized among all six SNP types, accounting for
Table 2
Summary of the Luciola leii transcriptome generated with Trinity assem-
bler software.
Item
Number
Table 1
Number of clean reads
Total transcripts
Total unigene length(nt)
Total unigene number
GC content (%)
294,391,102
105,644
54,231,639
54,282
41.24%
2409
Summary of RNA-Seq information for each sample.
Category
Control
C_1
BaP treatment
T_1 T_2
C_2
Total clean reads
75,866,242 71,280,880 74,922,326 72,321,654
Unigene N50 (nt)
Total mapped reads (%)
Unique match (%)
Multi-position match (%)
Total unmapped reads (%) 11.67
Perfect match (%)
Q20 percentage
88.33
64.90
23.43
88.51
63.99
23.62
12.39
48.36
97.53
83.46
67.47
20.99
11.54
47.81
97.47
83.61
65.16
22.37
12.47
48.07
97.51
Mean unigene length (nt)
NR annotation
Swiss-Prot annotation
GO annotation
KEGG annotation
All annotated unigenes
999
28,562
22,156
14,047
20,511
29,297
48.54
97.51

230
Q.-L. Zhang et al. / Science of the Total Environment 661 (2019) 226–234
A
2
2
1
1
5.0
0.0
5.0
0.0
2
1.3
1
8.9
17.7
1
0.9
1
0.4
6
.6
4
.4
5
0
.0
.0
1
.2
1.0 1.5 1.3
1.1
1.3
0
.1
0.5 0.4
0.5
0.1 0.1 0.1 0.4 0.1 0.1
B
2
2
1
1
50,000
00,000
50,000
00,000
5
0,000
0
Fig. 2. (A) Frequency of classified repeat SSR types. (B) SNP variation types of L. leii larvae.
high mapping rate of clean reads indicates the results of RNA sequenc-
ing were robust and can be used for further quantification of gene ex-
pression. Sample repeatability between two replicates was assessed
by calculating Pearson's correlation coefficient based on the FPKM
values of the unigenes. The Pearson's correlation coefficients were
N0.90 between the two biological replicates from treatment or control
samples (Fig. 1B and C). A high correlation between biological replicates
indicates that an effective experimental treatment was obtained and
that accidental error in a single sample was effectively avoided in this
study.
A total of 2414 DEGs were identified in response to BaP stress, of
which expression of 1350 DEGs were up-regulated, while 1064 DEGs
were down-regulated in the L. leii larvae treated by BaP in comparison
to control (Table S2). According to log2 FC in treatment compared to
the control, the top 15 up- and down-regulated DEGs are presented in
Table 3. The gene encoding mitochondrial arginyl tRNA synthetase
oxygen species (ROS) in fish liver, which significantly enhanced
hypoxia-related gene expression (Yu et al., 2008). Mitochondria, an en-
ergy production and metabolic factory, play the most important role in
hypoxia resistance of insects (Zhang et al., 2017c). In the list of the top
15 DEGs, several DEGs were detected that involved mitochondrial func-
tion, such as genes encoding superoxide-generating NADPH oxidase
heavy chain subunit C (NOXC), cytochrome c oxidase subunit 5A, mito-
chondrial (COX5A) and cytochrome-c1-heme lyase (CYT2). Desouki
et al. (2005) reported that NOX could generated ROS in human cells
and revealed that knockdown of mitochondrial genes caused down-
regulation of gene encoding NOXs while expression of NOX genes was
rescued by complementation with wild type mitochondrial genes. Pre-
vious studies reported that the oxidative stress caused by BaP disturbed
the expression of COX genes in earthworms (Eisenia fetida) (Lin et al.,
2008). Moreover, changes in expression of COX genes are sensitive,
and these genes have been demonstrated as an oxygen sensor in brain
neurons (Susann et al., 2010). Cytochromec1, which is a subunit of the
mitochondrial ubiquinol-cytochrome-c reductase, must be connected
with heme protein under CYT2 catalysis during cytochrome c1 import
into the mitochondria (Nicholson et al., 1989). Thus, we speculated
that ROS generation caused by BaP stress leads to differential expression
of these four function- and hypoxia-related mitochondrial genes (RARS,
NOXC, COX5A, CYT2) between BaP-treated and control groups. Further-
more, they may be potential biomarkers for oxidative stress following
xenobiotic contaminant exposure in aquatic L. leii. CYPs have been ex-
tensively studied in insects, with many investigations focusing on the
gene function of CYP in aquatic animals in response to BaP challenge.
(
RARS) showed the most up-regulated expression of the DEGs, followed
by cuticle protein 8 (CP8) and cytochrome P450 monooxygenase
CYP304E2 (CYP304E2). Conversely, the gene encoding cAMP-binding
protein 1 (CAPA1) was the most up-regulated gene, followed by the
gene encoding autophagy-related protein
8
(ATG8) and
chymotrypsin-1 (CHYM1). RARS is a gene that participates in protein
translation and hypoxia response. Previous research has demonstrated
that functional changes of RARS alters the survival rate of Caenorhabditis
elegans under hypoxic stress and that the tolerance level to hypoxia is
inversely correlated with the translation rate (Anderson et al., 2009).
Furthermore, BaP exposure induced rapid accumulation of reactive

Q.-L. Zhang et al. / Science of the Total Environment 661 (2019) 226–234
231
For example, Tian et al. (2014) demonstrated that mollusk CYPs can be
regulated by marine contaminants and the expression suppression of
several CYP-like genes may be caused by the high toxic effect and oxida-
tive damage of BaP in Chlamys farreri. In addition, studies show that BaP
suppression of immune system function relies on CYP-catalyzed pro-
cesses of immunotoxic metabolism in Japanese medaka (Oryzias latipes)
protein (HS718), heat shock protein 100 (HSP100) and DNA mismatch
repair protein Mlh1 (MLH1). This indicates that prevention of DNA
damage caused by BaP is a key biological process in the response of L.
leii to BaP stress. Additionally, we also identified two genes encoding cy-
tosolic carboxypeptidase 1 (AGTPBP1) and luciferin 4-monooxygenase
(LUCI). Roles of AGTPBP1 and LUCI are not well discussed due to com-
plexity and diversity of toxicological response, but these two genes
were identified as DEGs after BaP treatment in this study, thereby show-
ing the potential impact of BaP on their biological function. The genes
without annotated information are presented as hypothetical proteins,
and their function should be further analyzed.
(
Pannetier et al., 2019). Furthermore, CYPs have been widely consid-
ered as a set of suitable molecular indicators for BaP pollution
Camatini et al., 1998; Ryeo-Ok et al., 2013). In this study, expression
of genes encoding cytochrome P450 monooxygenase CYP304E2
CYP304E2), cytochrome P450 6BQ13 (CYP6BQ13) and cytochrome
(
(
P450 345B1 (CYP345B1) was altered by BaP exposure. Thus, our inves-
tigation suggests that BaP stress influences the detoxification function
of L. leii and also indicates that these CYPs may be useful as molecular
toxicological indicators of aquatic fireflies in water environments.
Interestingly, insect cuticle is the main component of the body wall
and functions by maintaining the morphology of the body, inhibiting
evaporation of body water and resisting invasion of xenobiotics
3.4. Functional enrichment analysis of differentially expressed genes
Following GO functional enrichment analysis of the DEGs, 46 biolog-
ical process terms were identified that involved the primarily immune
response, response to stimulus, nucleic acid metabolism and repair,
morphogenesis and development, and amino acid metabolic processes,
such as response to stimuli, immune response, DNA metabolism and the
cuticle chitin metabolic process (Table S3). For the eight enriched GO
terms associated with molecular function, we found that these catego-
ries were primarily related to biomacromolecule binding and structural
constituents of extracellular matrices and cuticles, such as sequence-
specific DNA binding, structural constituent of cuticle, extracellular ma-
trix structural constituent. Additionally, for the 13 cellular component
terms enriched by DEGs, those involving the extracellular region and
complex and nucleus function were overrepresented. In addition, 32
KEGG pathways enriched by the DEGs were obtained, indicating that
the gene expression of pathway members was altered under BaP stress
(Table S4) and that these pathways are key signaling transduces of a
progressive molecular response to BaP exposure in L. leii larvae. The sig-
nificantly enriched KEGG terms were primarily associated with xenobi-
otic and drug metabolism, disease, immune response, substance
digestion and absorption, such as protein digestion and absorption,
(Vincent and Wegst, 2004). Loss or reduction of cuticle protein content
leads to abnormal or even loss of cuticular function, thus affecting nor-
mal insect growth and development (Moussian et al., 2006). Therefore,
differential expression of the genes encoding cuticle protein 8 (CP8) and
cuticle protein 7 (CP7) suggest that BaP exposure probably affects the
developmental processes of L. leii larvae. Furthermore, we found up-
regulation of genes involved in ubiquitin and cell death, such as those
encoding ubiquitin-conjugating enzyme E2G (UBE2G), ubiquitin-like
protein 5 (UBL5) and cell death protein 2 (CED2) under BaP stress.
Song et al. (2019) found that various doses of BaP significantly affected
the expression of these related genes when using RNA sequencing to
study polar cod (Boreogadus saida). In addition, many studies have
found that heat shock proteins (HSPs) could protect against DNA dam-
age caused by BaP (Gong et al., 2007; Zheng et al., 2008). In this study,
several of the top 15 DEGs were related to the protection of DNA and re-
pair of DNA damage, such as hsp70/hsp90 organizing protein-like
Table 3
Top 15 up- and down-regulated genes between the control and treatment groups.
Gene ID
log2Ratio (treatment/control)
FDR
Annotation
Gene name
The top 15 most up-regulated genes
Unigene26775
Unigene10953
Unigene11072
Unigene10015
Unigene11983
Unigene33539
Unigene10425
Unigene8226
Unigene11616
Unigene8917
Unigene17694
Unigene19571
Unigene37979
Unigene29414
Unigene9415
17.04
15.87
15.25
15.18
14.93
14.72
14.59
14.14
13.81
13.71
13.71
13.50
13.48
13.46
13.44
7.37E-263
1.39E-223
2.81E-77
1.04E-107
3.36E-105
2.76E-56
1.27E-206
3.25E-52
4.85E-113
9.50E-102
2.47E-50
2.18E-24
2.17E-24
6.78E-22
6.01E-41
Mitochondrial arginyl tRNA synthetase
Cuticle protein 8
Hypothetical protein
Cytochrome P450 monooxygenase CYP304E2
Cytochrome P450 6BQ13
Ubiquitin-conjugating enzyme E2G
Hsp70/Hsp90 organizing protein-like protein
Hypothetical protein
Hypothetical protein
Cuticle protein 7
Cytochrome P450 345B1
Superoxide-generating NADPH oxidase heavy chain subunit C
Hypothetical protein
RARS
CP8
/
CYP304E2
CYP6BQ13
UBE2G
HS718
/
/
CP7
CYP345B1
NOXC
/
UBL5
CED2
Ubiquitin-like protein 5
Cell death protein 2
The top 15 most down-regulated genes
Unigene27586
Unigene30398
Unigene18974
Unigene12583
Unigene12275
Unigene31925
Unigene2849
Unigene20059
Unigene11399
Unigene24204
Unigene6022
Unigene6510
Unigene7788
Unigene36211
Unigene30489
−16.94
−14.22
−13.45
−13.08
−12.64
−12.44
−12.30
−12.15
−11.87
−11.78
−11.75
−11.54
−11.51
−11.35
−11.34
0
Hypothetical protein
/
3.33E-55
1.44E-22
8.76E-58
1.44E-22
6.73E-13
4.75E-12
5.34E-22
5.34E-22
1.14E-08
1.95E-06
8.56E-10
7.90E-08
1.32E-05
1.03E-06
cAMP-binding protein 1
Autophagy-related protein 8
Chymotrypsin-1
Cytochrome b5
Cytochrome c oxidase subunit 5A, mitochondrial
Hypothetical protein
Cytosolic carboxypeptidase 1
DNA mismatch repair protein Mlh1
Hypothetical protein
Luciferin 4-monooxygenase
Cytochrome c1, heme lyase, mitochondrial
Heat shock protein 100
CAPA1
ATG8
CHYM1
CYTB5
COX5A
/
AGTPBP1
MLH1
/
LUCI
CYT2
HSP100
/
Hypothetical protein
Hypothetical protein
/

232
Q.-L. Zhang et al. / Science of the Total Environment 661 (2019) 226–234
A
were also found in response to BaP stress in other animals, including
marine cephalochordate-amphioxus (Branchiostoma belcheri) (Zhang
et al., 2019a), coelenterate-coral (Acropora hyacinthus) embryo (Xiao
et al., 2018), teleostean polar cod (B. saida) (Song et al., 2019) and
annulata-earthworm (Eisenia fetida) (Zhang et al., 2017a). This analysis
indicates that, to a great extent, the aquatic firefly shares similar molec-
ular mechanisms regarding the BaP response with other animal groups.
ko00983
8
ko00982
1
0
8
3.5. Analysis of key BaP-responsive genes
0
4
Three well-documented KEGG terms involving xenobiotic biodegra-
1
1
dation and metabolism were significantly enriched, including metabo-
lism of xenobiotics by cytochrome P450, drug metabolism-
cytochrome P450 and drug metabolism-other enzymes (Guo et al.,
ko00980
2
013; Guo et al., 2016). A total of eight DEGs were shared among
B
these three pathways (Fig. 3A), suggesting a key role of these eight
DEGs in xenobiotic biodegradation. These genes were annotated to the
8
6
4
2
0
CYP
members,
glucuronosyltransferase
(UGT),
UDP-
glucuronosyltransferase 2B10-like (UGT2B10L), prostaglandin-H2 D-
isomerase (PTGDS), alcohol dehydrogenase, propanol-preferring
(
ADHP) and aldehyde dehydrogenase (NAD(P)+) (ALDH). These gene
family/genes encode several proteins that can metabolize a wide variety
of aquatic xenobiotic contaminants, such as nitrite, farm chemicals, or-
ganic compounds, heavy metal, plant-derived toxic substance and the
products of oxidative stress (i.e. free radicals and ROSs) (Guo et al.,
2
013; Guo et al., 2016; Zhang et al., 2019b). Thus, they are identified
-
4
-2
0
2
4
6
8
as key BaP-responsive genes (Gene set 1, GS1) in the current study. In
addition, to conduct an effective qRT-PCR analysis to validate the results
of RNA sequencing, four additional DEGs were randomly selected for
further analysis. Results of qRT-PCR analysis of the 12 DEGs showed a
close correlation (Pearson correlation coefficient = 0.912, p b 0.01) of
fold change in expression between the RNA sequencing and qRT-PCR
results, thereby supporting the accuracy of the RNA sequencing and
analysis in this experiment (Fig. 3B).
Expression changes of the GS1 genes in the BaP treatment in com-
parison to the control are presented in detail in Fig. 4 for multiple
time points. The results show that the expression levels of UGT,
CYP3A, CYP9 and CYP6AS5 were significantly induced at all time points
BaP-post-exposure (6, 12, 24 and 48). Meanwhile, the gene encoding
UGT2B10L showed significant down-regulation of expression at 6, 12
and 24 hpe. Furthermore, previous investigations reported that these
-
-
2
4
Fold change (RNA-seq)
Fig. 3. Venn diagram of DEGs belonging to three KEGG pathways involving xenobiotic
biodegradation and response. (D) Correlation between the relative fold changes in
expression from RNA-Seq and qRT-PCR analyses. Three technical replicates were
conducted for each of the two biological replicates.
drug metabolism - other enzymes, drug metabolism - cytochrome P450,
metabolism of xenobiotics by cytochrome P450, tyrosine metabolism,
Toll-like receptor signaling pathway, Vibrio cholerae infection and ma-
laria. Most of the signaling pathways identified in the current study
A
B
C
D
6
4
.0
.0
60
40
1.5
1.0
15
10
Control
Treatment
Control
Treatment
Control
Treatment
Control
Treatment
2
.0
0
20
0.5
5
0
0
1.5
1.0
0
6
h
12h 24h
UGT
48h
6h
12h 24h
CYP3A
48h
6h
12h 24h
UGT2B10L
48h
6h
12h 24h
CYP9
48h
E
F
TreatmentG
H
15
1
.5
Control
Treatment
3
0
0
Control
Treatment
Control
Control
Treatment
10
5
1
.0
2
0
.5
0
1
0
0.5
0
0
0
6
h
12h 24h
ALDH
48h
6
h
12h 24h
CYP6AS5
48h
6h
12h 24h
PTGDS
48h
6h
12h 24h
ADHP
48h
Fig. 4. qRT-PCR results of DEGs shared by the three xenobiotic biodegradation and metabolism pathways in L. leii larvae exposed to BaP stress at different time points. (A) UGT, (B) CYP3A,
⁎⁎
⁎
(
C) UGT2B10L, (D) CYP9, (E) CYP6AS5, (F) PTGDS, (G) ADHP and (H) ALDH. Data are presented as the mean ± SD (n = 6). p b 0.05, p b 0.01 relative to the corresponding time points of
the control; one-way ANOVA plus Bonferroni post-tests.

Q.-L. Zhang et al. / Science of the Total Environment 661 (2019) 226–234
233
Fu, X., Ballantyne, L., 2006. Luciola leii sp. nov., a new species of aquatic firefly (Coleoptera:
Lampyridae: Luciolinae) from mainland China. Can. Entomol. 138, 339–347.
Fu, X.H., Nobuyoshi, O., Lei, C.L., 2004. Morphological and biological observations on
aquatic firefly Luciola substriata (Gorham) (Coleoptera: Lampyridae) in China. Acta
Entomol. Sin. 47, 372–378.
Fu, X.H., Wang, Y.Y., Lei, C.L., 2005. Adaptive external morphology and swimming behav-
ior in the aquatic firefly. Luciola substriata. Entomol. Knowl. 42, 401–406.
Fu, X., Ballantyne, L.A., Lambkin, C.L., 2010. Aquatica gen. nov. from mainland China with a
description of Aquatica wuhana sp. nov. (Coleoptera: Lampyridae: Luciolinae).
Zootaxa 2530 (1), 18.
detoxifying proteins (CYPs and UGTs) participated in the response to
oxidative stress in rat astrocytes (Gradinaru et al., 2012). The gene
encoding alcohol dehydrogenase, propanol-preferring (ADHP) was sig-
nificantly up-regulated at 12, 24 and 48 hpe. A high expression of ADH
genes was identified as important to the growth of Metarhizium acridum
under hypoxia exposure by detoxification of acetaldehyde (Nora et al.,
2
011). In addition, the expression of PTGDS and ALDH genes were
down-regulated at different time points following BaP exposure. Sev-
eral studies have reported that the expression level of ALDH can be al-
tered by BaP stress, and ALDH is widely considered as a common
marker of cancer stem cells due to its significant effect on the expression
levels of stem cell markers in human BEAS-2B cells (Liu et al., 2016).
These findings indicate that oxidative stress caused by BaP exposure is
a crucial factor for triggering the metabolic and toxicological response
of L. leii larvae to xenobiotics in the water environment.
Gong, Z.Y., Niu, P.Y., Wang, F., Yuan, J., Zou, Y.F., Yang, J., Tang-Chun, W.U., 2007. Heat
shock protein 70 protects against DNA damage of A549 cells caused by benzo(a)
pyrene. J. Environ. Occup. Med. 24, 190–195.
Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Xian, A., Fan, L.,
Raychowdhury, R., Zeng, Q., 2011. Trinity: reconstructing a full-length transcriptome
without a genome from RNA-Seq data. Nat. Biotechnol. 29, 644–652.
Gradinaru, D., Minn, A.L., Artur, Y., Minn, A., Heydel, J.M., 2012. Effect of oxidative stress
on UDP-glucuronosyltransferases in rat astrocytes. Toxicol. Lett. 213, 316–324.
Guo, H., Ye, C.X., Wang, A.L., Xian, J.A., Liao, S.A., Miao, Y.T., Zhang, S.P., 2013. Trascriptome
analysis of the Pacific white shrimp Litopenaeus vannamei exposed to nitrite by RNA-
seq. Fish Shellfish Immunol. 35, 2008–2016.
Guo, H., Xian, J.A., Wang, A.L., 2016. Analysis of digital gene expression profiling in hemo-
cytes of white shrimp Litopenaeus vannamei under nitrite stress. Fish Shellfish
Immunol. 56, 1–11.
4
. Conclusions
Haas, B.J., Papanicolaou, A., Yassour, M., Grabherr, M., Blood, P.D., Bowden, J., Couger, M.B.,
Eccles, D., Li, B., Lieber, M., MacManes, M.D., Ott, M., Orvis, J., Pochet, N., Strozzi, F.,
Weeks, N., Westerman, R., William, T., Dewey, C.N., Henschel, R., LeDuc, R.D.,
Friedman, N., Regev, A., 2013. De novo transcript sequence reconstruction from
RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc.
For the first time, the toxicological response of aquatic insects to con-
taminants was explored using RNA sequencing. The transcriptome of
the aquatic firefly L. leii was generated based on large-scale sequencing
data. A number of molecular markers (SSRs and SNPs) that can be used
in investigations of genetic diversity were predicted, and DEGs involv-
ing BaP stress were also identified. These DEGs are primarily involved
in diverse signaling pathways related to disease, xenobiotic biodegrada-
tion, immune response and prevention of DNA damage. The expression
levels of eight key BaP-responsive genes, which were newly identified
in this study, were further investigated using qRT-PCR, thereby
confirming their crucial roles in response to BaP exposure in L. leii.
This study provides valuable information to further elucidate the molec-
ular response of aquatic insects following exposure to various environ-
mental stresses and indicates a complex and diverse network of
molecular toxicology in L. leii. However, due to the majority of unknown
factors related to the toxicological response to BaP exposure, several ob-
tained results that cannot be reasonably discussed require further
exploration.
8, 1494.
Jha, A.R., Zhou, D., Brown, C.D., Kreitman, M., Haddad, G.G., White, K.P., 2016. Shared ge-
netic signals of hypoxia adaptation in Drosophila and in high-altitude human popula-
tions. Mol. Biol. Evol. 33, 501–517.
Jia, X., Song, C., Chen, J., 2012. Advances in the study of toxicological effects of benzopy-
rene on fish immune system. Chin. Agric. Sci. Bull. 28, 110–115.
Langmead, B., Trapnell, C., Pop, M., Salzberg, S.L., 2009. Ultrafast and memory-efficient
alignment of short DNA sequences to the human genome. Genome Biol. 10, R25.
Li, X., Nair, A., Wang, S., Wang, L., 2015. Quality control of RNA-Seq experiments. Bioinfor-
matics 28, 2184–2185.
Lin, Z.S., Fang, S.Y., Yan, Q.X., Sun, T.H., Wei, Z., Ackland, M.L., 2008. Impact of benzo[a]
pyrene the expression of mitochondrion-encoded genes in the earthworm Eisenia
fetida. Environ. Sci. 29, 391–397.
Liu, Y., Lu, R., Gu, J., Chen, Y., Zhang, X., Lan, Z., Hao, W., Hua, W., Zeng, J., 2016. Aldehyde
dehydrogenase 1A1 up-regulates stem cell markers in benzo[a]pyrene-induced ma-
lignant transformation of BEAS-2B cells. Environ. Toxicol. Pharmacol. 45, 241–250.
Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using real-
time quantitative PCR and the 2−
ΔΔCT
method. Methods 25, 402–408.
McKenna, A., Hanna, M.E., Sivachenko, A., Cibulskis, K., Kernytsky, A., Garimella, K.,
Altshuler, D., Gabriel, S., Daly, M., DePristo, M.A., 2010. The genome analysis toolkit:
a MapReduce framework for analyzing next-generation DNA sequencing data. Ge-
nome Res. 20, 1297–1303.
Supplementary data to this article can be found online at https://doi.
org/10.1016/j.scitotenv.2019.01.156.
Miller, P.A., Elliott, N.G., Vaillancourt, R.E., Kube, P.D., Koutoulis, A., 2014. Genetic diversity
and pedigree assignment in tetraploid Pacific oysters (Crassostrea gigas). Aquaculture
Acknowledgments
433, 318–324.
Moussian, B., Seifarth, C., Müller, U., Berger, J., Schwarz, H., 2006. Cuticle differentiation
This study was supported by the Natural Science Foundation of
China (31760042 and 31760713).
during Drosophila embryogenesis. Arthropod Struct. Dev. 35, 137–152.
Nicholson, D.W., Stuart, R.A., Neupert, W., 1989. Biogenesis of cytochrome c1. Role of cy-
tochrome c1 heme lyase and of the two proteolytic processing steps during import
into mitochondria. J. Biol. Chem. 264, 10156–10168.
References
Nora, G., Srisombat, P., Macdonald, J.M., Gamcsik, M.P., Ngo, L.Y., Hohl, T.M., Cramer, R.A.,
2011. In vivo hypoxia and a fungal alcohol dehydrogenase influence the pathogenesis
Abdel-Shafy, H.I., Mansour, M.S.M., 2016. A review on polycyclic aromatic hydrocarbons:
source, environmental impact, effect on human health and remediation. Egypt. J. Pet.
of invasive pulmonary aspergillosis. PLoS Pathog. 7, e1002145.
Pannetier, P., Morin, B., Clérandeau, C., Lacroix, C., Cabon, J., Cachot, J., Danion, M., 2019.
Comparative biomarker responses in Japanese medaka (Oryzias latipes) exposed to
benzo[a]pyrene and challenged with betanodavirus at three different life stages.
Sci. Total Environ. 652, 964–976.
25, 107–123.
Anders, S., Huber, W., 2010. Differential expression analysis for sequence count data. Ge-
nome Biol. 11, R106.
Anderson, L.L., Xianrong, M., Scott, B.A., Michael, C., 2009. Survival from hypoxia in C.
elegans by inactivation of aminoacyl-tRNA synthetases. Science 323, 630–633.
Behringer, D., Zimmermann, H., Ziegenhagen, B., Liepelt, S., 2015. Differential gene ex-
pression reveals candidate genes for drought stress response in Abies alba (Pinaceae).
PLoS One 10, e0124564.
Borchert, J., Karbe, L., Westendorf, J., 1997. Uptake and metabolism of benzo(a)pyrene
absorbed to sediment by the freshwater invertebrate species Chironomus riparius
and Sphaerium corneum. Bull. Environ. Contam. Toxicol. 58, 158–165.
Camatini, M., Bonfanti, P., Colombo, A., Urani, C., 1998. Molecular approaches to evaluate
pollutants. Chemosphere 37, 2717–2738.
Pertea, G., Huang, X., Liang, F., Antonescu, V., Sultana, R., Karamycheva, S., Lee, Y., White, J.,
Cheung, F., Parvizi, B., Tsai, J., Quackenbush, J., 2003. TIGR Gene Indices clustering
tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformat-
ics 19, 651–652.
Ryeo-Ok, K., Bo-Mi, K., Chang-Bum, J., Nelson, D.R., Jae-Seong, L., Jae-Sung, R., 2013. Ex-
pression pattern of entire cytochrome P450 genes and response of defensomes in
the benzo[a]pyrene-exposed monogonont rotifer Brachionus koreanus. Environ. Sci.
Technol. 47, 13804–13812.
Song, Y., Nahrgang, J., Tollefsen, K.E., 2019. Transcriptomic analysis reveals dose-
dependent modes of action of benzo(a)pyrene in polar cod (Boreogadus saida). Sci.
Total Environ. 653, 176–189.
Susann, H., Cordian, B., Susanne, A., 2010. Effect of hypoxia on the transcription pattern of
subunit isoforms and the kinetics of cytochrome c oxidase in cortical astrocytes and
cerebellar neurons. J. Neurochem. 99, 937–951.
Conesa, A., Götz, S., García-Gómez, J.M., Terol, J., Talón, M., Robles, M., 2005. Blast2GO: a
universal tool for annotation, visualization and analysis in functional genomics re-
search. Bioinformatics 21, 3674–3676.
Daehwan, K., Ben, L., Salzberg, S.L., 2015. HISAT: a fast spliced aligner with low memory
requirements. Nat. Methods 12, 357–360.
Tang, P.A., Wu, H.J., Hao, X., Ju, X.R., Wei, S., Zhang, Q.L., Yuan, M.L., 2017. Characterization of
transcriptome in the Indian meal moth Plodia interpunctella (Lepidoptera: Pyralidae)
and gene expression analysis during developmental stages. Gene 622, 29–41.
Theissinger, K., Falckenhayn, C., Blande, D., Toljamo, A., Gutekunst, J., Makkonen, J., Jussila,
J., Lyko, F., Schrimpf, A., Schulz, R., 2016. De Novo assembly and annotation of the
freshwater crayfish Astacus astacus transcriptome. Mar. Genomics 28, 7–10.
Desouki, M.M., Kulawiec, M., Bansal, S., Das, G.C., Singh, K.K., 2005. Cross talk between mi-
tochondria and superoxide generating NADPH oxidase in breast and ovarian tumors.
Cancer Biol. Ther. 4, 1367–1373.
Fu, X.H., 2005. The Biology and Behavior of Two Species of Aquatic Fireflies in Mainland
China. Huazhong Agricultural University, Wuhan.

234
Q.-L. Zhang et al. / Science of the Total Environment 661 (2019) 226–234
Thiel, T., Michalek, W., Varshney, R., Graner, A., 2003. Exploiting EST databases for the de-
velopment and characterization of gene-derived SSR-markers in barley (Hordeum
vulgare L.). Theor. Appl. Genet. 106, 411–422.
Tian, S., Pan, L., Hui, Z., 2014. Identification of a CYP3A-like gene and CYPs mRNA expres-
sion modulation following exposure to benzo[a]pyrene in the bivalve mollusk Chla-
mys farreri. Mar. Environ. Res. 94, 7–15.
Triinu, K., Maido, R., 2007. Enhancements and modifications of primer design program
Primer3. Bioinformatics 23, 1289–1291.
Vincent, J.F.V., Wegst, U.G.K., 2004. Design and mechanical properties of insect cuticle. Ar-
Zhang, Q.L., Zhang, L., Yang, X.Z., Wang, X.-T., Li, X.-P., Wang, J., Chen, J.Y., Yuan, M.L.,
2017b. Comparative transcriptomic analysis of Tibetan Gynaephora to explore the ge-
netic basis of insect adaptation to divergent altitude environments. Sci. Rep. 7, 16972.
Zhang, Q.L., Zhang, L., Zhao, T.X., Wang, J., Zhu, Q.H., Chen, J.Y., Yuan, M.L., 2017c. Gene se-
quence variations and expression patterns of mitochondrial genes are associated
with the adaptive evolution of two Gynaephora species (Lepidoptera: Lymantriinae)
living in different high-elevation environments. Gene 610, 148–155.
Zhang, Q.L., Zhu, Q.H., Xie, Z.Q., Xu, B., Wang, X.Q., Chen, J.Y., 2017d. Genome-wide gene
expression analysis of amphioxus (Branchiostoma belcheri) following lipopolysaccha-
ride challenge using strand-specific RNA-seq. RNA Biol. 14, 1799–1809.
Zhang, Q.L., Qiu, H.Y., Liang, M.Z., Luo, B., Wang, X.Q., Chen, J.Y., 2017e. Exploring gene ex-
pression changes in the amphioxus gill after poly(I:C) challenge using digital expres-
sion profiling. Fish Shellfish Immunol. 70, 57–65.
thropod Struct. Dev. 33, 187–199.
Vongsangnak, W., Chumnanpuen, P., Sriboonlert, A., 2016. Transcriptome analysis reveals
candidate genes involved in luciferin metabolism in Luciola aquatilis (Coleoptera:
Lampyridae). PeerJ 4, e2534.
Wang, M.L., Barkley, N.A., Jenkins, T.M., 2009. Microsatellite markers in plants and insects.
Part I: applications of biotechnology. J. Genet. Genomics 29, 623–627.
Xiao, R., Zhou, H., Chen, C.M., Cheng, H., Li, H., Xie, J., Zhao, H., Han, Q., Diao, X., 2018. Tran-
scriptional responses of Acropora hyacinthus embryo under the benzo(a)pyrene
stress by deep sequencing. Chemosphere 206, 387–397.
Zhang, Q.L., Dong, Z.X., Xiong, Y., Li, H.W., Guo, J., Wang, F., Deng, X.Y., Chen, J.Y., Lin, L.B.,
2019a. Genome-wide transcriptional response of microRNAs to the benzo(a)pyrene
stress in amphioxus Branchiostoma belcheri. Chemosphere 218, 205–210.
Zhang, Q.L., Zhang, G.L., Xiong, Y., Li, H.W., Guo, J., Wang, F., Deng, X.Y., Chen, J.Y., Wang,
Y.-J., Lin, L.-B., 2019b. Genome-wide gene expression analysis reveals novel insights
into the response to nitrite stress in gills of Branchiostoma belcheri. Chemosphere
218, 609–615.
Yu, R.M.K., Ng, P.K.S., Tan, T., Chu, D.L.H., Wu, R.S.S., Kong, R.Y.C., 2008. Enhancement of
hypoxia-induced gene expression in fish liver by the aryl hydrocarbon receptor
(
AhR) ligand, benzo[a]pyrene (BaP). Aquat. Toxicol. 90, 235–242.
Zheng, S.L., Sun, T.H., Xiao, H., Qiu, X.Y., Song, Y.F., 2008. Low dose benzo(a)pyrene up-
regulated the transcription of HSP70 and HSP90 in Eisenia fetida. Chin. J. Appl. Ecol.
19, 401–406.
Zhang, L., Duan, X., He, N., Chen, X., Shi, J., Li, W., Xu, L., Li, H., 2017a. Exposure to lethal
levels of benzo[a]pyrene or cadmium trigger distinct protein expression patterns in
earthworms (Eisenia fetida). Sci. Total Environ. 595, 733–742.