Abstract
Epacromius coerulipes (Ivanov) is one of the most widely distributed
locusts. To date, the main methods to kill locusts still rely on
chemical controls, which can result in the selection of locusts with
resistance to chemical pesticides. Butene-fipronil is a new pesticide
that was discovered by the structural modification of fipronil. This
pesticide has been used to control various agricultural pests and has
become an important pesticide product to control pests that exhibit
resistance to other pesticides, including locusts. To extend its useful
half-life, studies of the initiation and progression of resistance to
this pesticide are needed. Herein, two E. coerulipes strains, a
pesticide-sensitive (PS) and a pesticide-resistant (PR) strain, were
chosen to undergo de novo assembly by paired-end transcriptome Illumina
sequencing. Overall, 63,033 unigenes were detected; the average gene
length was 772 bp and the N50 was 1,589 bp. Among these unigenes,
∼25,132 (39.87% of the total) could be identified as known proteins in
bioinformatic databases from national centers. A comparison of the PR
and PS strains revealed that 2,568 genes were differentially expressed,
including 1,646 and 922 genes that were up- and down-regulated,
respectively. According to the Gene Ontology (GO) database, among
biological processes the metabolic process group was the largest group
(6,900 genes, 22.47%) and contained a high frequency of differentially
expressed genes (544 genes, 27.54%). According to the Clusters of
Orthologous Groups (COG) categories, 28 genes, representing 2.98% of
all genes, belonged to the group of genes involved in the biosynthesis,
transportation, and catabolism of secondary metabolites. The
differentially expressed genes that we identified are involved in 50
metabolic pathways. Among these pathways, the metabolism pathway was
the most represented. After enrichment analysis of differential gene
expression pathways, six pathways—ribosome; starch, and sucrose
metabolism; ascorbate and aldarate metabolism; drug
metabolism-cytochrome P450; metabolism of xenobiotics by cytochrome
P450; and glutathione metabolism—showed a high degree of enrichment.
Among these pathways, drug metabolism-cytochrome P450, metabolism of
xenobiotics by cytochrome P450, and glutathione metabolism have been
associated with pesticide metabolism. Furthermore, 316 unigenes in the
E. coerulipes transcriptome encode detoxifying enzymes and 76 unigenes
encode target proteins of pesticides. Among these genes, 23 genes that
encode detoxifying enzymes in the resistance group were found to be
up-regulated. The transcriptome sequencing results of E. coerulipes
established a genomics database of E. coerulipes for the first time.
This study also establishes a molecular basis for gene function
analysis of E. coerulipes. Moreover, it provides a theoretical resource
for mechanistic studies on pesticide resistance through the screening
and investigation of resistance genes.
Keywords: pesticide resistance, Epacromius coerulipes (Ivanov),
transcriptome, differential gene expression analysis
__________________________________________________________________
Heilongjiang is one of the 10 provinces that contain grasslands in
China, and it contains 4.33 million hm^2 of grassland. However, ∼1.33
million hm^2 of grassland are attacked by locusts annually and
0.66–0.87 million hm^2 of grassland become disaster areas ([36]Zhou
et al. 2010). Recently, changes in the climate and cropland environment
in northern China have aggravated the regional occurrence of locusts
and resulted in an increasing trend for various locust together hazards
([37]Hong et al. 2014). Previous studies of the prairie locust in
Heilongjiang reported that the annual growth of dominant locust species
in Heilongjiang prairies, such as Epacromius coerulipes (Ivanov) and
Oxya chinensis Thunberg, can reach a population density of 300–500
insects/m^2 and greater. The occurrence rate of Epacromius coerulipes
accounts for up to 30% of the total, and exceeds 50% in some regions,
with serious and frequent occurrence each year and occasional outbursts
([38]Liu et al. 2008, [39]Zhou et al. 2011, Wang et al. 2014).
E. coerulipes, an insect belonging to Orthoptera, Acridoidea,
Edipodidae, is one of the locusts that is widely distributed in the
grasslands of China and in areas of both farming and herding. It not
only harms grasslands, as they have many hosts, and particularly thrive
on gramineous forage grass. They can grow in vegetation-sparse
wasteland and saline-alkali soil, degenerate grassland, and aggravate
salt and alkali conditions ([40]Tian et al. 2009); this can establish a
vicious circle of insect proliferation and damage to the ecological
balance in the prairie. E. coerulipes also damages gramineous crops and
broad-leaved crops, such as soybeans and alfalfa ([41]Tian et al. 2003,
[42]Tian et al. 2010, [43]Xu et al. 2013). Usually, hazard to crops
begins from the edge of a field, then expands to the middle and even
damages all crops in serious occurrences. Meanwhile, E. coerulipes can
fly short distances in cases of food shortage and cause widespread and
severe economic damage to agriculture and husbandry.
Because locust control in China focuses on chemical pesticides, locusts
have began to acquire chemical pesticide resistance and the efficiency
of locust control decreases throughout the continuous occurrence of
locust plagues. Consequently, it is important to study the molecular
mechanisms of pesticides on locusts to best use chemical pesticides and
extend their life span to sustainably control locusts. However, there
are few studies on pesticide resistance genes in locusts. Studies of
pesticide resistance, aided by transcriptome technology, have been
restricted to migratory locusts rather than E. coerulipes. Therefore,
high-flux Illumina sequencing techniques have been applied to discover
pesticide resistance genes by assembly, functional annotation, and
analysis, which has greatly improved our knowledge of E. coerulipes
pesticide resistance from a molecular perspective by studying
resistance mechanisms to provide a theoretical basis for mixing and
rotationally using chemical pesticides. In addition, such efforts may
allow the frequency of resistance genes to be limited as a consequence
of improved chemical prevention technology for locusts. Meanwhile, such
types of databases can provide useful data for the study of the genes
and genome of E. coerulipes.
Materials and Methods
Sampling and Feeding of Insects
Insects
E. coerulipes are collected from the natural prairie in Honggang
District, Daqing, Heilongjiang, with sensitivity to butene-fipronil of
LD[50] = 0.0303 μg/larva. They are subcultured in insect rearing cages
(length, height, and width of 1 by 1 by 1 m; cages were covered with
screens) under conditions of 28 ± 2°C, 65 ± 6% relative humidity, and a
photoperiod of 14:10 (L:D) h and fed fresh wheat seedlings.
Sensitive Strain Selection
The purpose was to develop sensitive strains by nonpesticide
treatments. After subculturing generations F5 and F7 in insect rearing
rooms, sensitivity to butene-fipronil reached LD[50] = 0.0298 μg/larva
and LD[50] = 0.0299 μg/larva, respectively. If an E. coerulipes strain
was highly sensitive to butene-fipronil, it is classified as a
sensitive strain for the following experiment, hereinafter referred as
PS.
Butene-Fipronil Resistance Strain Selection
If insects were sprayed with butene-fipronil and fed wheat seedlings,
their death rate ranged from 30 to 70% for resistance strain selection.
Compared with sensitive strains, generation F5 LD[50] = 0.1773 μg/larva
yielded a resistance ratio of 5.95, while the F7
LD[50] = 0.3499 μg/larva with a resistance ratio of 11.74 after
toxicity tests were conducted for subculture selection in insect
rearing rooms. In such cases, F7 will be applied as the resistance
strain for later experiments, hereinafter referred as PR.
Extraction and Testing of RNA
Fifth-instar larvae female nymphs of the above two strains (PS and PR
strain) were used for RNA extraction. Total RNA was extracted following
the manufacturer’s instructions (Invitrogen Trizol Reagent of RNA
extraction kit 15596018). This was carried out to ensure that the
samples for transcriptome sequencing were of a sufficient quality. The
purity, concentration, and integrity of RNA samples were tested
according to the manufacturer’s methods for the gel and Nanodrop 2000
(Thermo Fisher Scientific) and Agilent Bioanalyzer 2100 (Agilent).
cDNA Library Construction and Sequencing
When the quality of RNA samples was confirmed, using a NEB kit for
library preparation, a cDNA library was constructed as follows:
magnetic beads coupled with oligo(dT) were applied to enrich mRNA in
eukaryotes. The mRNA transcripts were randomly broken by bivalent
cations under high temperature conditions (94°C). The first cDNA chain
was synthesized using a target fragment of mRNA as a template and
randomly synthesized hexamers as primers. Then, buffer solution, dNTPs,
RNase H, and DNA polymerase I were added to synthesize the second cDNA
chain. The cDNA was next purified with AMPure XP beads. Purified
double-stranded cDNA was repaired at the end, and then an A-tail was
added and ligated with a sequencing joint. The fragment size (the
actual size was 150–250 bp without a joint) was selected using AMPure
XP beads (cDNA purification with 1.8 times the magnetic beads, i.e.,
100 µl reaction system and 180 µl magnetic beads). The cDNA library was
obtained after PCR enrichment (PCR conditions: 98°C 10 s, 65°C 30 s,
72°C 30 s, 72°C 5 min, 12 cycles). The prepared DNA library was
subjected to sequencing using an Illumina HiSeq2500 (Beijing Biomarker
Technologies Co.) with a read length of PE125. Two samples are
sequenced in the same RUN, but a different index was added to
distinguish them before sequencing.
The Transcriptome Assembly
Sample sequences were assembled to indirectly deepen the sequencing
results and obtain a more complete transcriptome dataset. To obtain
high-quality clean data, raw sequencing data were processed to remove
sequencing joints, poly-N regions, and low-quality reads. Meanwhile,
the Q20, Q30, GC content, and repetitive sequences were calculated and
all downstream analyses were carried out using high-quality clean data
(Q30 > 85%). High-quality sequencing data were assembled using Trinity
software ([44]Grabherr et al. 2011). The sequencing reads were first
broken into smaller fragments (K-mer), and then these fragments were
extended to generate a longer fragment (Contig). The fragment
collections (Component) were obtained based on overlaps of these
fragments. Transcript sequences were then recognized in each fragment
collection using a De Bruijin graph-based method and sequencing reads
data ([45]Grabherr et al. 2011, [46]Haas et al. 2013).
Functional Annotation of Genes
Unigene sequences were compared in the following databases using BLAST
software ([47]http://blast.ncbi.nlm.nih.gov/Blast.cgi) to obtain
annotation information for unigenes. Based on BLAST parameters, those
unigenes with an E-value less than 10^-5 were selected. Sequences of
selected unigenes were aligned within databases, including the NCBI
nonredundant protein sequences
(Nr: [48]ftp://ftp.ncbi.nih.gov/blast/db/), Clusters of Orthologous
Groups (COG: [49]http://www.ncbi.nlm.nih.gov/COG/), a manually
annotated and reviewed protein sequence database (Swiss-Prot:
[50]http://www.uniprot.org/), the Kyoto Encyclopedia of Genes and
Genomes (KEGG: [51]http://www.genome.jp/kegg/), and Gene Ontology (GO:
[52]http://www.geneontology.org/) databases.
Calculating the Expression Level of Unigenes
Sequencing reads of each sample were compared with the unigene database
using the Bowtie method ([53]Langmead et al. 2009). The comparison
results and RSEM ([54]Li et al. 2011) (RSEM is software used to
calculate levels) were used to estimate the expression levels with the
FPKM value (FPKM represents the expression level of genes) representing
the gene expression level.
[MATH: FPKM=cDNA
FragmentsMapped Fragments
(Millions)×Transcript Length
(kb) :MATH]
In the above formula, cDNA fragments represent the number of aligned
fragments within certain transcripts and indicate the number of double
reads. Mapped fragments (Millions) represent the total number of
aligned fragments within transcripts with a unit of 10^6. Transcript
length (kb) is the length of transcripts with a unit of 10^3 bases.
Then, information was obtained for each gene.
Detection and Cluster Analysis of Differentially Expressed Genes
Differentially expressed gene sets for the two samples were acquired by
differential expression analysis using EBSeq ([55]Leng et al. 2013). In
the differential gene expression analysis, significant p-values based
on the original hypothesis was corrected using the Benjamini–Hochberg
method ([56]Anders et al. 2010). The corrected p-value, i.e., the False
Discovery Rate (FDR), was used as a key index to screen differentially
expressed genes to minimize the false-positive results caused by
independent statistical hypothesis tests for high-value gene expression
values. During this screening process, FDR < 0.01 and Fold Change
(FC) ≥ 2 were used as screening criteria. FC represented the ratio of
expression values between two samples. The selected differentially
expressed genes were analyzed for hierarchical clustering ([57]Murtagh
et al. 2014). Identically or similarly expressed genes were clustered
to reveal the various expression patterns under different experimental
conditions.
Enrichment Analysis of Differentially Expressed Genes
After be identified by screening,the differentially expressed genes
were analyzed by GO function enrichment, COG classification, pathway
annotation analysis, and KEGG pathway enrichment analysis.
The differentially expressed genes identified by screening were used to
do GO function enrichment and draw histogram: extract GO annotation of
differential expression gene, map GO function to the corresponding
secondary features on the background of all unigene’s GO annotation
([58]Ye et al. 2006), and then obtain accurate results of GO function
enrichment by using the Fisher’s exact test.
COG classification was used to extract the Cog annotation results of
differential expression gene, and then drawing histogram to select
([59]http://www.ncbi.nlm.nih.gov/COG/).
KEGG pathway enrichment analysis by KOBAS2.0
([60]http://kobas.cbi.pku.edu.cn/home.do) ([61]Xie et al. 2011)
extracted the pathway information of differential expression gene, and
then on the background of all unigene’s pathway information, calculated
the enrichment factor and detected significance by the Fisher’s exact
test. When enrichment factors were smaller and the absolute Q value was
higher, the enrichment level of a differentially expressed gene was
more significant and reliable. The calculation formula for the
Enrichment Factor was as follows:
[MATH: Enrichment
Factor=Total number of differentially expressed genes/Total gene number in
the KEGGNumber of differentially expressed
genes in a pathway/Total gene number in a
pathway :MATH]
Results
Sequencing Analysis and Assembly
The extracted cDNA samples from E. coerulipes were subjected to
sequencing using an Illumina sequencing platform. The data shown in
[62]Table 1 indicate that the sequencing results were highly accurate
and could be utilized for future studies. The sequencing assembly was
performed using Trinity software ([63]Grabherr et al. 2011), and the
parameters are shown in [64]Table 2. These findings suggest that the
sequencing assembly was ideal.
Table 1.
Summary of the raw reads for the two samples
Samples Read number Clean Data GC content Cycle Q20% %≥Q30
PR (T01) 26,869,615 6,766,002,460 50.71% 100.00 90.63%
PS (T02) 26,428,337 6,653,605,161 50.81% 100.00 90.09%
[65]Open in a new tab
Note: GC content: Clean Data G and C percentage of the total bases;
Cycle Q20: Average quality score is greater than or equal to 20% of the
cycle, Q30: Quality Score of base is greater than or equal to 30% of
the total bases.
Table 2.
Summary of Illumina transcriptome assembly for E. coerulipes
Length range Contig Transcript Unigene
200–300 6,232,095 (99.35%) 30,249 (31.46%) 26,218 (41.59%)
300–500 15,971 (0.25%) 18,617 (19.36%) 14,150 (22.45%)
500–1,000 10,968 (0.17%) 16,364 (17.02%) 9,772 (15.50%)
1,000–2,000 7,402 (0.12%) 14,974 (15.57%) 6,855 (10.88%)
2,000+ 6,195 (0.10%) 15,947 (16.59%) 6,038 (9.58%)
Total number 6,272,631 96,151 63,033
Total length 306,171,970 105,999,072 48,663,063
N50 length 47 2,269 1,589
Mean length 48.81 1,102.42 772.03
[66]Open in a new tab
Note: The sequencing assembly obtained 6,272,631 Contigs and the N50
length was 47.96, 151 transcripts were obtained. The total length was
105.999 Mb, the average length was 1,102 bp, and the N50 length was
2,269. Among them, the number of transcripts with a length over 1 kb
was 30,921, which was 32.16% of the total number. Additionally, 15,947
of transcripts with a length over 2 kb were identified, which was
16.59% of all transcripts. After clustering and assembly analysis of
the transcripts, 63,033 unigenes were obtained, with a total length of
48.663 Mb and an average length of 772 bp, and a N50 length of
1,589 bp. Among them, 12,893 unigenes with a length over 1 kb were
identified, which represents 20.45% of the total, and 6,083 unigenes
with a length over 2 kb were identified, which represents 9.58% of the
total.
Gene Functional Annotation
As shown in [67]Table 3, 25,132 unigenes were successfully annotated,
which represents 39.87% of all genes. Among these genes, there were
24,841 Nr database annotated unigenes, which represents 39.41% of all
genes. The Swissprot database contained 16,490 annotated unigenes,
which represents 26.16% of total genes. The COG database had 8,031
annotated unigenes, which represents 12.71% of all genes. The GO
database included 11,558 annotated unigenes, which represents 18.34% of
total genes. The KEGG database had 7,218 annotated unigenes, which
represents 11.46% of all genes. However, 37,901 unigenes, which
represents 60.13% of all genes, were unannotated. These findings
indicate that these unigenes can yield short reads by sequencing
technology, and can be subjected to various types of analysis ([68]Hou
et al. 2011).
Table 3.
Functional annotation of the E. coerulipes transcriptome
Annotated database Annotated unigenes number 300 ≤ length < 1,000
Length ≥ 1,000
COG_Annotation 8,013 (12.71) 2,603 3,955
GO_Annotation 11,558 (18.34%) 3,863 5,064
KEGG_Annotation 7,218 (11.46%) 2,260 3,694
Swissprot_Annotation 16,490 (26.16%) 5,676 7,889
nr_Annotation 24,841 (39.41%) 9,079 9,881
All_Annotated 25,132 (39.87%) 9,185 9,895
[69]Open in a new tab
Note: Nr database annotated 24,841 unigenes, which represents 39.41% of
the total. The Swissprot database annotated 16,490 unigenes, which
represents 26.16% of the total. The COG database annotated 8,031
unigenes, which represents 12.71% of the total. The GO database
annotated 11,558 unigenes, which represents 18.34% of the total. The
KEGG database annotated 7,218 unigenes, which represents 11.46% of the
total.
Annotation and Clustering Analysis of Differentially Expressed Genes
Statistical analyses and annotations for differentially expressed genes
between the PS and PR strains of E. coerulipes are presented in
[70]Table 4, among which total 2,568 (4.07%) differentially expressed
genes were acquired, and proportion of total annotations in each
database reached 7.15%. Up-regulated genes accounted for 64.10% of all
differentially expressed genes. The blue and red zones, respectively,
represent differences with low and high levels of significance (Fig.
1).
Table 4.
Number of differentially expressed genes and gene annotation
DEGs number Annotated COG GO KEGG Swiss-Prot nr
All DEG 2,568 1,798 760 815 584 1,424 1,784
Up-regulated 1,646 1,177 537 579 453 956 1,168
Down-regulated 922 621 223 236 131 468 616
[71]Open in a new tab
Note: The 2,568 differentially expressed genes were obtained and 1,798
of them were annotated to each database. Up-regulated genes included
1,646 and 1,177 genes from those that were annotated; down-regulated
genes included 922 and 621 genes from those that were annotated.
Functional Enrichment Analysis of Differentially Expressed Genes
The GO database contains three main categories of annotation: molecular
function, cellular component, and biological process ([72]Fig. 2). In
the cellular component group, the 3,000 unigenes in a cell or cellular
organelle were the largest category. The number of differentially
expressed genes was also the greatest, and both the number of unigenes
in a cell or cellular organelle were greater than 270. In the molecular
function group, the total number of catalytic activity and
differentially expressed genes from them were greatest and represented
the highest proportion of all genes. The proportion of differentially
expressed genes increased in certain subgroups, including structural
molecular activity, electron carrier activity, antioxidant activity,
and channel regulator activity. In the biological processes group, the
number of genes involved in metabolic processes was the greatest
(6,870, representing 27.34% of all genes). The number of differentially
expressed genes was also the greatest and had the highest proportion
(544, 30.25%) of all genes. The proportion of differentially expressed
genes (132, 0.43%; 10, 0.51%) also increased in the immune system
processes subgroup.
Fig. 1.
[73]Fig. 1.
[74]Open in a new tab
Clustering gene expression pattern of differentially expressed genes.
Note: Each column represents one sample, each row represents one gene,
and colors represent the expression level, FPKM with log[2] in sample
genes. The blue area represents a low significant difference and the
red area represents a highly significant difference.
Fig. 2.
[75]Fig. 2.
[76]Open in a new tab
GO classification. Annotation statistics of differentially expressed
genes in the secondary node of GO. Note: The horizontal axis shows
secondary nodes of three categories in GO. The vertical axis displays
the percentage of annotated genes versus the total gene number. The red
columns display annotation information of the total genes and the blue
columns represent annotation information of the differentially
expressed genes.
In [77]Fig. 3, general function prediction alone was the largest class
group; the number of differentially expressed genes was 133, which
represents 14.16% of all genes. The number of differentially expressed
genes involved in carbohydrate transport and metabolism represents
12.89%; in translation, ribosomal structure, and biogenesis represents
12.14%; the number of differentially expressed genes in other subgroups
represents <10%. Because secondary metabolites in insects significantly
affect insecticides, secondary metabolite biosynthesis, transport, and
catabolism is an important class group among the COG categories. The
number of differentially expressed genes was 28, which represents 2.98%
of the total genes.
Fig. 3.
[78]Fig. 3.
[79]Open in a new tab
COG function classification of consensus sequences. Note: The
categories of the COG are shown on the horizontal axis and gene numbers
are plotted on the vertical axis. The number of differentially
expressed genes involved in general function prediction is 133; in
carbohydrate transport and metabolism is 121; in translation, ribosomal
structure, and biogenesis is 114; in amino acid transport and
metabolism is 84; in posttranslational modification, protein turnover,
and chaperones is 79; in replication, recombination, and repair is 61;
in inorganic ion transport and metabolism is 56; in energy production
and conversion is 5; in lipid transport and metabolism is 47; and in
cell wall/membrane/envelope biogenesis is 30; the number of
differentially expressed genes in other subgroups is <30.
As shown in [80]Fig. 4, the KEGG annotation results of differentially
expressed genes were classified according to KEGG pathway
classification. The differentially expressed genes were distributed
across 50 metabolic pathways. Among them, the metabolism pathway
included the greatest number of genes. The number of differentially
expressed genes involved in the metabolism of xenobiotics by cytochrome
P450, drug metabolism of cytochrome P450, drug metabolism of other
enzymes, and glutathione metabolism, were 14, 15, 7, and 15,
respectively. These metabolic pathways are all related to the
degradation of pesticides in insects.
Fig. 4.
[81]Fig. 4.
[82]Open in a new tab
KEGG categories of differentially expressed genes. Note: The vertical
axis lists the names of the metabolic pathways in the KEGG, and the
horizontal axis shows the proportion of annotated genes in each pathway
versus the total number of annotated genes. The differentially
expressed genes are distributed in 50 metabolic pathways, the pathway
of metabolism holds the highest all number of genes.
KEGG enrichment analysis of differentially expressed genes is shown in
[83]Fig. 5, and 20 selected metabolic pathways were analyzed. This
analysis indicated that six metabolic pathways are of the highest
concentration when diagrams of the top 20 metabolic pathways with high
concentration were analyzed. These six pathways—ribosome, starch and
sucrose metabolism; ascorbate and aldarate metabolism; drug
metabolism-cytochrome P450; metabolism of xenobiotics by cytochrome
P450; and glutathione metabolism—showed the highest enrichment scores.
Among these pathways, three pathways—drug metabolism-cytochrome P450;
metabolism of xenobiotics by cytochrome P450; and glutathione
metabolism—are related to pesticide metabolism.
Fig. 5.
[84]Fig. 5.
[85]Open in a new tab
Scatterplot of differentially expressed genes from the KEGG pathway
enrichment analysis. Note: Every graph in the figure represents one
KEGG pathway and the name of the pathway is shown on the right side of
the legend. The horizontal axis shows the enrichment factor, which
represents the proportion of the annotated differentially expressed
genes in each pathway compared to the total annotated gene number in
that pathway. The smaller enrichment factor indicates a higher
enrichment score of differentially expressed genes in the pathway. The
vertical axis displays an absolute Q value, which is the p-value tested
by multiple hypothesis and calibrations. Thus, a larger value on the
vertical axis indicates a more reliable enrichment score of
differentially expressed genes in a pathway.
Annotation and Analysis of Pesticide Resistance Genes
In the E. coerulipis transcriptome, many gene sequences that encoded
pesticide detoxification enzymes and target proteins were annotated. In
this present study, 316 genes that encoded pesticide detoxification
enzymes were annotated, including 43 glutathione transferase (GST)
genes, 90 carboxylesterase (CarEs) genes, and 183 cytochrome P450
genes. There were 74 unigenes that encoded target proteins of
pesticides, including 10 γ-aminobutyric acid receptor (GABA), 15
nicotinic acetylcholine receptor (nAChRs), 3 ryanodine receptor, 38
AChE enzyme, and 10 voltage-gated sodium channel (VGSC) genes.
Additionally, 23 genes, including 5 GST, 8 CarEs, and 10 cytochrome
P450s, were up-regulated in the PR strain compared with the PS strain
([86]Table 5).
Table 5.
Unique transcripts associated with insecticide resistance in E.
coerulipes
Enzymes or target Annotated unigenes number Up-regulated number
(PS-vs.-PR)
GSTs 43 5
CarEs 90 8
CYP450 183 10
GABA 10 –
AchR 15 –
RyR 3 –
AchE 38 –
CGSC 10 –
[87]Open in a new tab
Note: GSTs, glutathione S-transferases; CarEs, carboxylesterase;
CYP450, cytochrome P450; GABA, γ-aminobutyric acid receptor; nAchR,
nicotinic acetylcholine receptor; RyR, ryanodine receptor; AchE,
acetylcholinesterase; VGSC, voltage-gated sodium channel; Up-regulated
number (PS-vs.-PR), GSTs: c8747.graph_c0, c34411.graph_c0,
c33067.graph_c0, c40846. graph_c0, c31380.graph_c0; CarEs :
c33871.graph_c0, c30386.graph_c0, c33452.graph_c0, c22727.graph_c0,
c27991.graph_c0, c34612.graph_c0, c35381.graph_c0, c30559.graph_c0;
CYP450: c23816.graph_c0, c37409. graph_c0, c29735.graph_c0,
c34866.graph_c0, c34963.graph_c0, c34617. graph_c0, c39389.graph_c0,
c37409.graph_c0, c22350.graph_c0, c39389. graph_c0.
Discussion
As a consequence of increasing research efforts and the completion of
insect genomes and transcriptomes, the breadth and depth of molecular
knowledge on insects has significantly expanded. However, locust
genetics received relatively less attention. [88]Badisco et al. (2011a)
studied the central nervous system (CNS) of the desert locust
(Schistocerca gregaria gregaria) based on expressed sequence tags (EST)
reads and obtained 34,672 original EST sequences from the CNS. Those
ESTs were assembled into 12,709 transcript sequences; one-third of them
were annotated. Based on the EST database of S. gregaria gregaria,
[89]Badisco et al. (2011b) compared the CNS genes between solitary and
gregarious desert locusts, and they identified 214 differentially
expressed genes. [90]Chen et al. (2010) studied the transcriptomes of
the migratory locust (Locusta migratoria) by de novo assembly. This
study yielded a deeper understanding of the genetics of incomplete
metamorphosis in insects and the origins of their deformation.
Furthermore, understanding the genes and pathways involved in the
development of locusts and the phase transitions in locusts could be
helpful for controlling these insects. [91]Jiang et al. (2012) first
analyzed retroelements in the migratory locust transcriptome and
identified 105 retroelements; this finding allowed us to better
understand the retrotransposon distribution in the migratory locust
transcriptome. Moreover, these findings showed that non-long terminal
repeat (LTR) retrotransposons in migratory locust transcriptomes were
extremely rich and exhibited great diversity. [92]Lv (2012) and
[93]Yang (2013) analyzed the transcriptomes of nymphs and adults of O.
chinensis and those of three insects strains of Atractomorpha sinensis,
respectively. These studies enhanced the transcriptome database of
orthopteran insects. Transcriptome sequencing was used to investigate
the resistance in insects, such as L. migratoria, the silverleaf
whitefly (Bemisia tabaci), the brown planthopper (Nilaparvata lugens),
and the greenhouse whitefly (Trialeurodes vaporariorum) ([94]Wang 2009,
[95]Wang et al. 2010, [96]Xue et al. 2010, Karatolos et al. 2012).
These studies greatly enhanced the efficiency and quantity of gene
annotation ([97]Ansorge. 2009) and improved research into insect
resistance to insecticides. These findings initially resulted in the
application of the transcriptome technique to study E. coerulipes
resistance to butane-fipronil. After butane-fipronil resistance
selection, E. coerulipes can acquire 2,568 (1,178 are annotated)
differentially expressed genes, which are distributed across 50
metabolic pathways; the largest category of these pathways was
metabolic pathways. Additionally, it has been reported that the drug
metabolism-cytochrome P450, metabolism of xenobiotics by cytochrome
P450, and glutathione metabolism pathways are abundantly expressed and
are closely linked to pesticide metabolism by insects, which provides a
theoretical basis for tracking metabolic resistance in detoxifying
enzymes.
Pesticide detoxification enzymes in insects, i.e., cytochrome P450s,
CarEs, and GSTs, participate in the metabolism of allogenic materials,
plant secondary substances, and pesticides ([98]Strode et al. 2008,
[99]Wang 2011). A total of 316 unigenes among the E. coerulipes
transcriptomes encode detoxifying enzymes. Different types and
quantities of pesticide detoxification enzyme genes have been reported
in previous studies based on the analysis of other insect
transcriptomes, such as 35 P450s, 6 CarEs, and 7 GST genes in B.
tabaci; 53 P450s, 49 CarEs, and 15 GST genes in Dialeurodes citri; and
102 P450s genes in Cimex hemipterus ([100]Wang 2011, [101]Mamidala
et al. 2012, [102]Chen et al. 2013). Furthermore, based on an analysis
of insect genome sequencing results, 85, 106, 164, and 83 P450s genes
have been found in Drosophila melanogaster, Anopheles gambiae, Aedes
aegypti, and Acyrthosiphon pisum (Harris), respectively ([103]Adams
et al. 2000, [104]Holt et al. 2002, [105]Strode et al. 2008,
[106]Ramsey et al. 2010).
Many studies have established that pesticide resistance in insects can
be achieved by the overexpression of various types of detoxification
enzymes. For example, the activity of P450 mono-oxygenase was greatly
increased in fipronil-resistant Sogatella furcifera ([107]Jian
et al.2010). Fipronil resistance in Bemisia tabaci and Laodelphax
striatellus is derived from increasing activity of p450 mono-oxygenase
([108]Kang et al. 2006, [109]Wang et al. 2010). Expression changes in
the carboxylesterase genes pxae22 and pxae31 can influence the
sensitivity of Plutella xylostella to fipronil ([110]Ren et al. 2015).
The fipronil resistance of Chilo suppressalis is related to the
activity and expression levels of carboxylesterase and p450
mono-oxygenase ([111]Li et al. 2007). Because butene-fipronil and
fipronil are the same kind of insecticide, they may have the same
resistance mechanism. The above analysis indicates that 23 newly
discovered up-regulated detoxifying enzymes may contribute to the
metabolism of butane-fipronil and show certain connections with the
onset of resistance. Consequently, it is necessary to verify those
genes and their expression levels to identify the exact relationship
between E. coerulipes resistance to butene-fipronil and detoxifying
enzymes in order to control resistance development by detoxifying
enzyme inhibitor for the efficient use of chemical pesticides to
achieve locust control.
The molecular targets of pesticides include the GABA receptor,
nicotinic acetylcholine receptor, ryanodine receptor, AChE enzyme, and
voltage-gated sodium channel (VGSC). These target proteins have been
reported to be associated with pesticide resistance. It has also been
reported that mutations in these genes encoding these target proteins
lead to varying degrees of insensitivity to pesticides in insects.
Mutations at M918V, L925I, and T929V in the VGSC gene have been
associated with pyrethroid resistance ([112]Chung et al. 2011).
Single-site mutations in a conserved position, Y151S, in two nicotinic
acetylcholine receptor subunit genes, Nla1 and Nla3, have been linked
to imidacloprid resistance in the brown planthopper ([113]Liu et al.
2005). The four mutated sites in the gene encoding the ryanodine
receptor have been associated with chlorantraniliprole resistance in
Plutella xylostella ([114]Guo et al. 2014). Endosulfan resistance in B.
tabaci ([115]Houndete et al. 2010) has been associated with mutations
in the GABA receptor subunit genes. The GABA receptor is also the
molecular target of phenyl pyrazole pesticides (butene-fipronil and
fipronil). Rdl gene mutation (Rdl1a→Rdl1s) of the P. xylostellas GABA
receptor is linked to their strong resistance to fipronil. The
fipronil-resistant Laodelphax striatellus ([116]Chen. 2012) may be
related to A2` mutations N in the GABA receptor gene. Our study found
that 76 unigenes encode target proteins of pesticides, and the
annotation of genes that encode target proteins in this present study
may represent the foundation for next-generation screening of mutation
sites.
Transcriptome research on E. coerulipes has rapidly and effectively
screened out resistance-related detoxifying enzymes and target site
genes, in addition to providing valuable data for further study of the
E. coerulipes resistance mechanism and control strategies.
Acknowledgments