Open Access

Comparative transcriptome analysis of lufenuron-resistant and susceptible strains of Spodoptera frugiperda (Lepidoptera: Noctuidae)

BMC Genomics201516:985

https://doi.org/10.1186/s12864-015-2183-z

Received: 4 July 2015

Accepted: 30 October 2015

Published: 21 November 2015

Abstract

Background

The evolution of insecticide resistance in Spodoptera frugiperda (Lepidoptera: Noctuidae) has resulted in large economic losses and disturbances to the environment and agroecosystems. Resistance to lufenuron, a chitin biosynthesis inhibitor insecticide, was recently documented in Brazilian populations of S. frugiperda. Thus, we utilized large-scale cDNA sequencing (RNA-Seq analysis) to compare the pattern of gene expression between lufenuron-resistant (LUF-R) and susceptible (LUF-S) S. larvae in an attempt to identify the molecular basis behind the resistance mechanism(s) of S. frugiperda to this insecticide.

Results

A transcriptome was assembled using approximately 19.6 million 100 bp-long single-end reads, which generated 18,506 transcripts with a N50 of 996 bp. A search against the NCBI non-redundant database generated 51.1 % (9,457) functionally annotated transcripts. A large portion of the alignments were homologous to insects, with the majority (45 %) being similar to sequences of Bombyx mori (Lepidoptera: Bombycidae). Moreover, 10 % of the alignments were similar to sequences of various species of Spodoptera (Lepidoptera: Noctuidae), with 3 % of them being similar to sequences of S. frugiperda. A comparative analysis of the gene expression between LUF-R and LUF-S S. frugiperda larvae identified 940 differentially expressed transcripts (p ≤ 0.05, t-test; fold change ≥ 4). Six of them were associated with cuticle metabolism. Of those, four were overexpressed in LUF-R larvae. The machinery involved with the detoxification process was represented by 35 differentially expressed transcripts; 24 of them belonging to P450 monooxygenases, four to glutathione-S-transferases, six to carboxylases and one to sulfotransferases. RNA-Seq analysis was validated for a number of selected candidate transcripts by using quantitative real time PCR (qPCR).

Conclusions

The gene expression profile of LUF-R larvae of S. frugiperda differs from LUF-S larvae. In general, gene expression is much higher in resistant larvae when compared to the susceptible ones, particularly for those genes involved with pathways for xenobiotic detoxification, mainly represented by P450 monooxygenases transcripts. Our data indicate that enzymes involved with the detoxification process, and mostly the P450, are one of the resistance mechanisms employed by the LUF-R S. frugiperda larvae against lufenuron.

Background

The fall armyworm Spodoptera frugiperda (Lepidoptera: Noctuidae) occurs mainly in tropical and subtropical regions [1], causing large losses to cotton and corn crops in the Americas [2]. Recent changes in the cropping system in the Brazilian savanah (Cerrado) by the integration of crops such as cotton, corn, soybean and millet have led to an increase in the population densities of S. frugiperda [3]. Control of S. frugiperda has been based on the use of genetically modified plants that express Bt (Bacillus thuringiensis) toxins and on the use of synthetic insecticides. However, the strong selective pressure caused by these control methods has led to an increase of S. frugiperda resistance to Bt toxins and to insecticides from different chemical groups [3, 4].

Chitin biosynthesis inhibitors act by interfering with the synthesis or deposition of chitin on the exoskeleton and on other chitinized structures of insects [5]. Given the specificity of their mode of action, this group of insecticides has been shown to have great potential in integrated pest management (IPM) programs because of its low toxicity to humans and higher animals [6]. The high pressure resulting from the widespread adoption of chitin biosynthesis inhibitors such as lufenuron for the control of insect pests in the Brazilian Cerrado has modified the susceptibility of S. frugiperda populations [7]. The observed reduction in the susceptibility of natural populations of S. frugiperda to chitin-synthesis inhibitors may indicate the evolution of resistance to this insecticide. Thus, measures that enable the preservation of the useful life of this molecule are necessary.

Currently, the main resistance mechanisms of S. frugiperda to insecticides involve mutations that lead to the insensitivity of target sites and/or alterations in the activity of enzymes involved with the detoxification or sequestration of xenobiotics [810]. Despite these findings, little is known about the gene expression profile of strains resistant to chitin biosynthesis inhibitors, and there are no studies for resistant strains of S. frugiperda to such products.

Next-generation sequencers (NGS), such as Solexa/Illumina™ (Illumina©), 454 (Roche©) and SOLID™ (Applied Biosystems©) platforms, generate a large amount of data [11] that allows for an in depth investigation of transcriptomes at a low cost and in a short period of time.

We characterized the larval transcriptome of S. frugiperda and compared the larval gene expression patterns between LUF-R and LUF-S strains as a step towards understanding the molecular mechanism(s) involved in S. frugiperda resistance to lufenuron in order to support the further development of rational, sustainable tools for pest and resistance management strategies.

Results

De novo assembly of a reference transcriptome

Sequencing of cDNA libraries in the Illumina HiScan 1000® platform generated 68,027,513 reads of approximately 100 bp, corresponding to 6,802,571,300 nucleotides. After selection of reads through quality filters, 23.2 % of the reads obtained were discarded due to their low quality scores (see Additional file 1).

Changes in the k-mer parameter led to changes in almost all observed variables (see Additional file 2). Thus, assemblies with a k-mer of 23, 25, 47, 53 or 55 were selected for performing the reference de novo assembly. Assemblies used 14,337,437 nucleotides, equivalent to 71.8 % of the bases submitted for analysis. These individual assemblies were merged in a single de novo reference assembly, which resulted in the generation of 18,506 transcripts with sizes ranging from 100 to 6,517 bp (see Additional file 3), with a mean length of 774.75 bp (see Additional file 4), a N50 of 996 bp and a N90 of 411 bp.

Functional annotation

From the 18,506 transcripts obtained, only 51.1 % (9,457) were functionally annotated after a heuristic search against the NCBI non-redundant protein database (see Additional file 5). The similarity analysis yielded e-values from 10−3 to 10−32 for nearly 30 % of the mapped transcripts, while 24 % transcripts had even lower e-values, ranging from 10−32 to 10−61. The highest e-value significance scores were identified in 6 % of the transcripts, and ranged from 10−148 to 10−177 (see Additional file 5). Almost all of the alignments obtained were related to insects (see Additional file 6), and their majority (45 %) aligned against the silk moth Bombyx mori (Lepidoptera: Bombycidae), followed by the monarch butterfly Danaus plexippus (Lepidoptera: Nymphalidae) (20 %). The representatives of the genus Spodoptera (Lepidoptera: Noctuidae) observed in the analysis were S. littoralis, S. exigua, S. litura and S. frugiperda, corresponding to 10 % of the obtained alignments; S. frugiperda accounted for 3 % of them. Nearly two-thirds (6,078) of the 9,457 annotated transcripts and gene ontologies were distributed in 13 functional groups in the category of proteins that function as cellular components, 11 groups in molecular functions and 23 groups in biological functions (see Additional file 7).

Differential gene expression between LUF-S and LUF-R S. frugiperda larvae

A comparative gene expression analysis demonstrated LUF-R and LUF-S S. frugiperda larvae have distinct gene expression profiles (Fig. 1). Among the 18,506 transcripts evaluated, 940 were differentially expressed (t-test, p ≤ 0.05; fold change ≥ 4) (Fig. 2). These transcripts are from a large variety of functional categories, with a great number of transcripts associated with catalytic, binding and metabolic processes (Fig. 3). Among the differentially expressed transcripts, 25 transcripts were overexpressed in the larvae of the LUF-R strain, with fold change values greater than 100 times the expression of the larvae of the LUF-S strain. The fold change of the top 20 highly expressed transcripts in the LUF-R larvae ranged from 125 to 620 fold the expression in the LUF-S larvae. Conversely, eight transcripts were suppressed in the LUF-R strain, varying from −110 to −1,974 fold the expression in the LUF-S strain (see Additional file 8).
Fig. 1

Distribution of differentially expressed transcripts from an RNA-seq analysis of susceptible (LUF-S) and lufenuron-resistant (LUF-R) strains of S. frugiperda, induced or non-induced by lufenuron

Fig. 2

Distribution of the transcripts of the S. frugiperda transcriptome based on the comparative analysis of the gene expression of susceptible (LUF-S) and lufenuron-resistant (LUF-R) strains, induced or non-induced by lufenuron. Marked in green are the transcripts with a significant difference in the expression level, based on the discriminative significance values (t-test, p < 0.05) and relative expression (>4) adopted here

Fig. 3

Distribution of gene ontology (GO) attributed to differentially expressed transcripts in susceptible (LUF-S) and lufenuron-resistant (LUF-R) S. frugiperda strains

The grouping analysis showed that the majority of the transcripts (61.3 %) were overexpressed, whereas 38.7 % were suppressed in the resistant strain when compared to the susceptible strain (Fig. 4). Our analysis did not indicate substantial differences in the pattern of expression within each strain (LUF-S and LUF-R) following their exposure to lufenuron, showing that the majority of the overexpressed and suppressed transcripts of the LUF-R strain exhibit constitutive expression.
Fig. 4

Comparative distribution of functionally annotated transcripts of lufenuron-susceptible (LUF-S) and lufenuron-resistant (LUF-R) S. frugiperda strains that showed changes in the expression level. RPKM values were represented as a scale of colors ranging from green to red, which will encompass values from the lowest (green) to the highest (red) RPKM values

Among the differentially expressed genes, two groups of candidate transcripts were selected because of their association with the cuticle metabolism and detoxification processes, as they could be related to the resistance mechanism of S. frugiperda to lufenuron (Figs. 5 and 6). Two out of the six selected transcripts involved in cuticle metabolism (L_663_T_4/6 and L721_T_4/5) were suppressed, while the remaining was overexpressed in the LUF-R as compared to LUF-S larvae (Fig. 5).
Fig. 5

Comparative distribution of transcripts associated with cuticle metabolism of lufenuron-susceptible (LUF-S) and lufenuron-resistant (LUF-R) S. frugiperda strains that showed changes in the expression level. RPKM values were represented as a scale of colors ranging from green to red, which will encompass values from the lowest (green) to the highest (red) RPKM values

Fig. 6

Comparative distribution of transcripts associated with detoxification enzymes of lufenuron-susceptible (LUF-S) and lufenuron-resistant (LUF-R) S. frugiperda strains that showed changes in the expression level. RPKM values were represented as a scale of colors ranging from green to red, which will encompass values from the lowest (green) to the highest (red) RPKM values. The title of each transcript consists of the identification of genes and code of the transcript. In P450’s transcribed capital letter and number is the class ID and group associated with gene superfamily

Out of the 35 differentially expressed genes involved in detoxification processes, 24 of them were annotated as P450 (CYP) monooxygenase, four as glutathione-S-transferase (GST), six as carboxylase (CCE) and one as sulfotransferase (SUR). CYP was represented by several family members, especially CYP9A, CYP3, CYP4 and CYP6. Transcript levels of 26 of these genes were much higher in LUF-R larvae, while nine of them had their expression drastically reduced (Fig. 6).

qPCR analysis

The differential expression patterns between LUF-S and LUF-R larvae were validated by the relative expression obtained by qPCR using a selected group of transcripts (Fig. 7). All of the transcripts had higher gene expression in the LUF-R as compared to the LUF-S larvae, especially the transcripts CYP9A9 - L_464_T_3/3 and CYP321A1- L_669_T_9/12, which exhibited fold changes of 45 and 900 times, respectively.
Fig. 7

qPCR analysis of selected CYP transcripts identified as differentially expressed in a broad RNA-Seq analysis of the larval transcriptome of susceptible (LUF-S) and resistant (LUF-R) strains of S. frugiperda to lufenuron, exposed (induced) or not (non-induced) to lufenuron treatment. Expression of the selected genes is provided as fold change (ΔΔCt) using their expression at the LUF-S, non-induced larvae as a reference

Discussion

Insect resistance is an evolutionary phenomenon arising from the continuous exposure of a population to the selective pressure represented by the indiscriminate use of insecticides [12]. Currently, the main mechanisms associated with resistance development are mutations in the target sites of insecticides, the intermediation of the metabolism of insecticides by detoxification enzymes and tegumental changes that limit insecticide penetration [13].

Next-generation sequencing (NGS) technologies have brought great advances for genomic studies in non-model organisms. These technologies provide a large amount of data at a low cost [14], increasing the possibility of recovering important biological information from transcriptomes [15]. Therefore, this information is of great importance, given that events such as insect resistance to insecticides are biologically complex phenomena related to adaptive processes, such as mutations and metabolic processes, that are vital for organism maintenance [12].

The de novo assembly of the transcriptome we performed allowed the evaluation of the differential expression between LUF-R and LUF-S larvae of S. frugiperda. A large number of differentially expressed transcripts were observed in between LUF-R and LUF-S larvae. To date, however, most of these transcripts have not been associated with the molecular mechanisms involved with insect resistance to insecticides. These results are similar to those reported for chlorantraniliprole-susceptible and resistant strains of Plutella xylostella (Lepidoptera: Plutellidae) [16]. In this case, 1,215 differentially expressed transcripts were identified, many of which were associated with metabolic pathways different from those clearly involved with the detoxification processes or of target sites of the insecticide. As an example, we observed severe down regulation of a eukaryotic translation initiation factor (L_194_T_3/3 transcript) in LUF-R as compared to LUF-S larvae, with a fold change of 1,974. Eukaryotic translation initiation factors have conserved functions related to translation events, mRNA recruitment and regulation of the cellular machinery for protein synthesis [17]. Therefore, significant changes in the expression of this transcript can be related to the strong regulation of post-transcriptional events in protein synthesis. Conversely, the L_1097_T_4/4 transcript, which is associated with the ubiquinol-cytochrome c reductase complex and linked to electron transport in the cellular respiration process [18], was the up-regulated transcript with the highest fold change (620) observed in the LUF-R strain. Regulation of genes that affect the post-transcriptional gene expression and post-translational protein modification, as the translational initiation factors and genes involved in protein ubiquitination, suggests that such mechanisms may also be involved in the observed resistance of S. frugiperda to lufenuron as demonstrated in the response of other organisms exposed to toxicants [1924].

Despite a great portion of the differentially expressed transcripts was not functionally characterized, our results reveal the complexity of adaptive processes resulting from the selection pressure from the continuous exposure to lufenuron. The increase in the number of genomic-based studies associated with insect resistance has strengthened the notion that insect resistance to various insecticides can be related to polygenic and/or epigenetic factors [25], as is the case of the strain we have analyzed [26].

There have been a number of studies looking at the mode of action and on the possible resistance mechanisms of insects to benzoylphenylureas, including lufenuron, but no mechanism has been clearly demonstrated yet [27]. Many studies have argued that resistance could arise from an elevated activity of the enzymes involved in chitin processing [5, 2729], but we did not observe significant changes in the expression of enzymes involved either in chitin synthesis (chitin synthases) or chitin degradation (chitinases). Some have also argued that resistance could be involved with the production of cuticle proteins [27]. In fact, we did detect up-regulation of transcripts associated with RR-1 cuticle proteins, L_1615_T_6/8 (cuticular protein isoform a), L_13_T_1/14 (cuticular protein rr-1 motif 46), L_13_T_11/14 (cuticular protein rr-1 motif 46) and L_623_T_2/4 (cuticle protein), but all of them are proteins associated with more flexible regions of the cuticle [30]. These results agree with those observed for diflubenzuron-treated Tribolium castaneum (Coleoptera: Tenebrionidae), as no changes in the expression levels of chitin synthases and chitinases were observed, although significant changes in the expression of cuticle proteins could be detected [27]. Our findings and the existing data available suggest that genes involved in chitin metabolism, modification and degradation are not the targets of benzoylphenylureas. Despite of the resistance of Tetranychus urticae to etoxazole was associated with mutations in conserved regions of the chitin synthase gene [31], we did not detect mutations in this gene in the LUF-R strain of S. frugiperda.

Alterations in the expression levels or mutations in the lufenuron receptor could also explain the resistance to lufenuron observed in the LUF-R strain of S. frugiperda. The lufenuron receptor remains to be described, but Abo-Elghar et al. [28] indicated the ABCC transporter sulfonylurea receptor (Sur) as the receptor for benzoylphenylurea due to the structural similarities of benzoylphenylurea and sulfonylurea. However, we did not detect any ABCC transporter or sulfonylurea receptor-like transcript in our transcriptomic analysis. Yet, the role of Sur as a sulfonylurea has been recently challenged, as this receptor was demonstrated to be dispensable for chitin synthesis in Drosophila melanogaster (Diptera: Drosophilidae), suggesting on the existence of an alternative sulfonylurea-sensitive ABC transporter to be involved with chitin synthesis and cuticle formation [32]. Therefore, we discarded this receptor as a possible mechanism of resistance in the LUF-R S. frugiperda strain.

In vitro studies have indicated benzoylphenylureas act through inhibition of the incorporation of N-acetylglucosamine (GlcNAc) into chitin in an ecdysteroid-dependent manner [33]. We did not find consistent changes in the expression of ecdysteroidogenic genes in between LUF-S and LUF-R strains of S. frugiperda that would suggest regulation of ecdysteroid synthesis, but we did observed a lower expression of ecdysteroid-22-kinase. This enzyme is involved in ecdysteroid inactivation in silkworm adult ovaries [34], suggesting that LUF-R larvae would retain higher titers of active ecdysteroids. Nevertheless, we did not observe any indication on differential production or incorporation of UDP-N-acetylglucosamine into chitin, as no changes in gene expression of UDP-N-acetylglucosamine diphosphorylases (GlcNAc production) or of chitin synthase (GlcNAc incorporation) were detected. Altogether, we did not find enough support to suggest the ecdysteroid-mediated regulation of GlcNAc synthesis or incorporation into chitin as a possible mechanism of resistance for the lufenuron-resistant S. frugiperda strain.

It has been demonstrated that benzoylphenylureas may affect cuticle formation in a concentration-dependent manner, with lower doses affecting cuticle thickness and microfibril orientation, while higher doses completely disrupt chitin synthesis [35]. The reported effects of benzoylphenylurea on the orientation and synthesis of chitin in the cuticle suggest an effect through disruption of the circadian clock involved in chitin secretion and deposition during cuticle formation [36], but we did not find evidence on the participation of clock genes as a mechanism of resistance to lufenuron in the LUF-R strain of S. frugiperda as well. We did though observe UDP-glucosyltransferases that were highly expressed in the LUF-R strain (L_1390_T_6/9 and L_1390_T_6/9), an enzyme that is involved with the process of chitin synthesis. However, recent progress due to large scale genomic sequencing has demonstrated UDP-glucosyltransferases form a multigenic family in insects [37, 38]. The diversity of UDP-glucosyltransferases in lepidopterans have also indicated their contribution in the process of detoxification through glycosylation [38, 39], and there are several indications in which an increase in the expression level of UGTs have been related to insect resistance to insects such DDT [40] and carbamate [41]. However, one other UDP-glucosyltransferase (L_492_T_3/8) highly similar to UGT40R3 was down regulated in the LUF-R strain, but UGT40R3 has been reported to be an UDP-glucosyltransferases exclusively associated with chemoreceptors present in insect antennae [42].

Our differential expression analysis between LUF-S and LUF-R S. frugiperda larvae identified a large number of transcripts associated with P450, GSTs and CCEs, of which many were overexpressed in LUF-R larvae. Detoxification of insecticides has been widely reported as one of the main mechanisms of insect adaptations to the high selective pressure exerted by insecticides [10]. Enzymes coded by genes from cytochrome P450, GSTs and CCEs are widely associated with the resistance mechanisms of insects to pesticides because of the degradation, detoxification and/or sequestration of xenobiotics [12]. The up-regulation of P450 enzymes and carboxylases, in addition to glycosyltransferases, sulfotransferases and glutathione-S-transferases, was recently demonstrated in a functional genomic analysis of a population of Tribolium castaneum contaminated with diflubenzuron, corroborating their association with the detoxification of xenobiotics [27].

A higher number of P450 monooxygenases was observed among the ESTs with higher relative expression. P450 monooxygenases have been one of the main classes of enzymes associated with lepidopteran resistance to insecticides such as pyrethroids [43], organophosphates [44] and diamides [16]. Several studies have been conducted aiming to evaluate the activity of these enzymes at the toxicological level [45]. Their role in insecticide degradation has been demonstrated by observed increased mortality of resistant populations of S. frugiperda to pyrethroids, organophosphates and carbamates when exposed mainly to monooxygenase inhibitors [46]. The role of CYPs in lufenuron detoxification has also been suggested, as CYP12A4 was up-regulated in a population of Drosophila melanogaster resistant to lufenuron, suggesting that the enzymes from cytochrome P450 play a key role in the resistance against this insecticide.

The finding of multiple overexpressed CYP families in LUF-R larvae also points for their role in the resistance of S. frugiperda to lufenuron. CYP3, CYP4, CYP6 and CYP9 are P450 families that have been argued as one of the mechanisms involved in lepidopteran resistance to insecticides [4749]. The high expression levels associated with the genes encoding detoxification enzymes, even in the absence of insecticide, show that the high expression of these genes occurs constitutively.

The metabolism of insecticides in insects certainly involves a series of complex metabolic processes. Even though the biochemical processes related to the detoxification are well described, there are important gaps related to the physiological and molecular mechanisms that govern the detoxification process [50]. We were able to demonstrate the overexpression of genes involved in the functionalization (phase I – eg CYPs) and conjugation (phase II – eg UDP-glucosyltransferases, GSTs), but not of genes involved in the excretion (phase III – eg ABC transporters) of xenobiotics in the process of lufenuron detoxification in the LUF-R strain of S. frugiperda [50]. This may indicate that the process of detoxification of lufenuron can mainly occur within the lumen of the gut.

The diverse gene expression profile between LUF-S and LUF-R strains indicate resistance also affects several other pathways, which may be related with the resistance fitness costs or may represent an additional mechanism to contribute with resistance of the LUF-R S. frugiperda to lufenuron. Therefore, there is the need for studies that provide a better understanding of the metabolic pathways adjacent to the classical detoxification pathways [12]. For many years, studies of the resistance mechanisms of insects to insecticides were scarce. However, the advance of sequencing technologies and the ability to generate large volumes of information have provided a better understanding of the mechanisms involved in the processes that lead to the resistance of insects to insecticides [16].

Conclusions

Gene expression of LUF-R larvae of S. frugiperda is much higher than of LUF-S in general, and particularly of those genes involved the detoxification process, including CYP, CCE and GST. We concluded that the high abundance of highly expressed CYP genes in the resistant as compared to the susceptible strain is one of the major mechanisms involved in resistance of S. frugiperda to lufenuron.

Methods

Insect preparation

Two laboratory-selected strains of S. frugiperda were used for sequencing, a lufenuron-susceptible strain (LUF-S) and a lufenuron-resistant strain (LUF-R). LC50s for LUF-S and LUF-R were estimated at 0.23 and 216.6 μg.mL−1 of lufenuron, respectively [26]. Larvae of fourth instars were fasted for 24 h. After this period, susceptible and resistant larvae were split into subgroups that were submitted or not to induction with lufenuron. Larvae were transferred to containers with an artificial diet treated superficially with 3.2 μg.mL−1 of lufenuron (Match 50 EC, 500 g.L−1 of lufenuron, emulsifiable concentrate, Syngenta Proteção de Cultivos Ltda.) diluted in distilled water added with 0.1 % Triton® and remained in contact with the food source for one hour (induced). In the case of control treatments, the artificial diet was surface-treated only with the 0.1 % Triton in water (non-induced). Therefore, larvae selected for RNA extraction and cDNA sequencing comprised four groups: induced and control LUF-S, and induced and control LUF-R.

Tissue sampling and total RNA extraction

Total RNA was extracted from five fourth instars per treatment using Trizol® Reagent (Invitrogen™), according to the manufacturer’s instructions. Samples were macerated in liquid nitrogen and homogenized in 1 mL Trizol reagent by vortexing, followed by the addition of 200 μL. Sample was mixed by vortexing and centrifuged for 15 min at 12,000 g. The aqueous phase was recovered, added to 500 μL of isopropanol and centrifuged (10 min x 12,000 g). The pellet obtained was washed in 1 mL of 75 % ice-cold ethanol, and the sample was centrifuged once again (5 min x 7,500 g). The pellet was air-dried at room temperature and resuspended in 60 μL of sterile milliQ water. Residual DNA was eliminated after treatment with 9 U of RNAse-free DNAse. The supernatant was recovered and added to 300 μL of Trizol and 200 μL of chloroform. Sample was agitated and centrifuged for 10 min at 13,000 g. The supernatant was discarded, the pellet was recovered in 300 μL of chloroform, and the solution was centrifuged for 5 min at 13,000 g. The supernatant was recovered, 1 mL of 100 % ethanol was added and the sample centrifuged again (10 min x 13,000 g). The pellet obtained was air-dried at room temperature and resuspended in DEPC-treated water. RNA quantity and quality were assessed in a Picodrop-Spectrophotometer 4.0.4.0.

The mRNA population in each sample was enriched by removing ribosomal RNA (rRNA) using the commercial system RiboMinus® Eukaryote Kit (Invitrogen™), following the manufacturer’s instructions. mRNA-enriched samples of LUF-S and LUF-R larvae of S. frugiperda were used for cDNA library preparation using the TruSeq Stranded mRNA Library Prep Kit, and sequencing in a HiSeq 1000® (Illumina©) platform at the Laboratory of Animal Biotechnology, Department of Animal Science, ESALQ, University of São Paulo, a sequencing service provider.

De novo assembly

The paired-end protocol was used and originated reads of approximately 100 bp for the sequencing of cDNA libraries. The reads obtained from the Illumina HiScan 1000® platform sequencing were filtered for removing bases with quality scores lower than 30 (where a Phred score of 30 corresponds to a 0.1 % expected error rate) at both the 5′and 3′ends and used for the assembly of a single reference transcriptome.

Then, aiming to maximize the computational performance, the duplicate reads were excluded, minimizing the time used for the de novo assembly of the transcriptome. The paired-end obtained were subjected to exploratory analysis using the Velvetoptimiser version 2.2.5 (https://github.com/tseemann/VelvetOptimiser) [51] to evaluate the diversity of transcripts assembled in different motif lengths (k-mers) [52]. Therefore, assemblies were obtained for all odd sizes of k-mers, varying from 19 to 61, following [51]. Assemblies with k-mers of 23, 25, 47, 53 and 55 obtained with the software Velvet version 1.2.10 (https://github.com/dzerbino/velvet/tree/master) were used to obtain a reference transcriptome. The concatenation of contigs of the assemblies obtained from the different k-mers used was performed in Oases version 0.2.08 (www.ebi.ac.uk/~zerbino/oases/) [53]. Transcripts sharing 95 % similarity were grouped using the CD-hit version 4.6 (https://github.com/weizhongli/cdhit) [54] to reduce the redundancy of the final assembly.

Functional annotation

The obtained transcripts were annotated after a similarity search against the NCBI non-redundant database using the BLASTx algorithm [55] available with the BLAST2GO software [56], with a cutoff value of 10−3. Enzyme classification (EC) codes and the annotation of metabolic pathways (KEGG - Kyoto Gene and Genomes) were generated from the direct mapping of the GO terms with their equivalent enzyme codes [57]. For functional classification, the consensus sequences were compared to those from the Interpro protein signature databases using the InterproScan with a cutoff of 10−5 [58].

Differential expression between the LUF-S and LUF-R S. frugiperda strains

The differential expression between the LUF-R and LUF-S strains of S. frugiperda was evaluated by determining the number of reads per kilobase of transcript per million reads (RPKM) mapped against the de novo assembled reference transcriptome earlier described, using only reads that counted against one target reference transcript following the remaining default parameters in the CLC Genomics Workbench software (QIAGEN Company). The experimental design comprised of two groups, the LUF-R or the LUF-S strain, each one counting with the subgroups induced and non-induced. To obtain a high specificity for the results, only the reads similar to a single transcript were computed, with a minimum similarity of 80 % and a read length alignment higher than 90 %. The mean RPKM values obtained for the resistant strain (induced and non-induced) were compared to those of the susceptible strain (induced and non-induced) by the t-test (p < 0.05), using the tools available in the CLC Genomics Workbench software (QIAGEN Company). Only transcripts with significant differences in the expression levels based on the discriminative significance values (p ≤ 0.05) and relative expression (fold change >4) were considered differentially expressed.

The validation of the comparative analysis using RNA-seq data between LUF-R and LUF-S S. frugiperda strains was performed for a set of differentially expressed transcripts by using quantitative real-time PCR (qPCR). Transcripts possibly associated with the mechanisms of resistance to lufenuron, such as those belonging to cuticle metabolism and metabolic detoxification enzymes were randomly selected as candidate genes (Table 1).
Table 1

Primer sets used in qPCR reactions for validation of the RNA-Seq analysis

Transcript

 

Sequence (5′- 3′)

L_464_T_3/3

Forward

CAATGCAATTCCTTGGACCAA

Reverse

GCACACCATCGATCCAATGA

L_406_T_12/23

Forward

TCGGTGAGAGGTATGCCAAAT

Reverse

AAGTTTCGCAAGACGTGCACTA

L_669_T_9/12

Forward

CAAACCAGCCTGCACCTGTA

Reverse

GGGCAACAGGACGTGTATAGG

L_1141_T_4/4

Forward

CCATCGAGGTTGTGCAAAAA

Reverse

CCTCAGCCAGTGTCCTCCAT

L_1141_T_1/4

Forward

CAATCCGCTCGCGTACATG

Reverse

AATCTTGACCCAATGCAATTCC

GAPDH

Forward

CGGTGTCTTCACAACCACAG

Reverse

TTGACACCAACGACGAACAT

qPCR validation of candidate genes

For validation of our differential expression analysis between LUF-S and LUF-R strains of S. frugiperda, 2 μg of RNA treated with RNAse-free DNAse obtained from each treatment was used in a reverse transcription reaction to produce cDNA. For the reverse transcription reaction, the commercial system ImProm-II™ Reverse Transcription System (Promega©) was used according to the manufacturer’s instructions. The cDNA was stored at −20 °C until use.

qPCR reactions consisted of 0.4 μg of cDNA, 12.5 μL of Maxima SYBR Green/ROX qPCR Master Mix® (Thermo Scientific©), 0.6 μM of each of the primers (Table 1) and 10.3 μL of water, totaling a final volume of 25 μL. The primers for each selected target were designed using the Primer Express® software (Life Technologies©), based on the sequences of the S. frugiperda transcriptome. The amplification reactions were performed in a ViiA™ 7 Real Time PCR System (Applied Biosystems®, Life Technologies©) platform. The thermocycling conditions used were as follows: 2 min of preheating at 50 °C, 10 min for initial denaturation at 95 °C, followed by 35 cycles at 95 °C for 15 s, 60 °C for 30 s and 72 °C for 30 s. The dissociation curves were analyzed to assess the specificity of the amplification. The standardization of the amplification was performed using glyceraldehyde 3-phosphate dehydrogenase (GAPDH) [59].

The experiment followed a completely randomized design, composed of four treatments similar to those used for the sequencing, containing three biological replicates (each replicate consisted of five larvae); each replicate was analyzed in technical triplicates. Differences in gene expression between samples analyzed were determined following the ΔΔCt method [60].

Availability of supporting data

All Illumina data have been deposited in NCBI’s Sequence Read Archive (SRA) under accession number PRJNA299878.

Abbreviations

Bt: 

Bacillus thuringiensis

CCE: 

carboxylase

CYP: 

cytochrome P450

DDT: 

1,1,1-trichloro-2,2-di(4-chlorophenyl)ethane

GAPDH: 

Glyceraldehyde 3-phosphate dehydrogenase

GO: 

Gene Ontology

GST: 

Glutathione S-transferase

IPM: 

Integrated Pest Management

KEGG: 

Kyoto Gene and Genomes

LC50: 

lethal concentration 50

LUF-R: 

Strain Resistance to lufenuron

LUF-S: 

Strain Susceptible to lufenuron

NCBI: 

National Center for Biotechnology Information

NGS: 

Next Generation Sequencing

qPCR: 

Quantitative real-time PCR

RPKM: 

Reads per kilo base per million mapped reads

SUR: 

Sulfonylurea receptor

UGT: 

UDP-glycosyltransferase

Declarations

Acknowledgements

This work was part of research conducted in partial fulfillment of the M.Sc. Degree in Entomology by ARBN at Escola Superior de Agricultura “Luiz de Queiroz” (ESALQ/USP). We thank National Council for the Improvement of Higher Education (CAPES) for the M.Sc. scholarship to ARBN and postdoctoral fellowship to PF. We also thank National Council for Scientific and Technological Development (CNPq) for granting research fellowships to FLC (Process 312094/2013-2) and CO (Process 312086/2013-0).

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Departamento de Entomologia e Acarologia, Escola Superior de Agricultura “Luiz de Queiroz”, Universidade de São Paulo

References

  1. Rojas JC, Virgen A, Malo EA. Seasonal and nocturnal flight activity of Spodoptera frugiperda males (Lepidoptera : Noctuidae) monitored by pheromone traps in the coast of Chiapas, Mexico. Fla Entomologist. 2004;87(4):496–503.View ArticleGoogle Scholar
  2. Meagher RL, Nagoshi RN. Population dynamics and occurrence of Spodoptera frugiperda host strains in southern Florida. Ecol Entomol. 2004;29(5):614–20.View ArticleGoogle Scholar
  3. Diez-Rodriguez G, Omoto C. Herança da resistência de Spodoptera frugiperda (J.E. Smith) (Lepidoptera: NOctuidae) a lambda-cialotrina. Neotrop Entomol. 2001;30:311–16.View ArticleGoogle Scholar
  4. Storer NP, Kubiszak ME, King JE, Thompson GD, Santos AC. Status of resistance to Bt maize in Spodoptera frugiperda: lessons from Puerto Rico. J Invertebr Pathol. 2012;110(3):294–300.View ArticlePubMedGoogle Scholar
  5. Merzendorfer H, Zimoch L. Chitin metabolism in insects: structure, function and regulation of chitin synthases and chitinases. J Exp Biol. 2003;206(24):4393–412.View ArticlePubMedGoogle Scholar
  6. Beeman RW. Recent advances in mode of action of insecticides. Annu Rev Entomol. 1982;27:253–81.View ArticlePubMedGoogle Scholar
  7. Schmidt FB. Baseline susceptibility of Spodoptera frugiperda (Lepidoptera: Noctuidae) to Lufenuron in corn. Piracicaba: Luiz de Queiroz College of Agriculture, University of São Paulo; 2002.Google Scholar
  8. Yu SJ, Nguyen SN, Abo-Elghar GE. Biochemical characteristics of insecticide resistance in the fall armyworm, Spodoptera frugiperda (J.E. Smith). Pestic Biochem Physiol. 2003;77(1):1–11.View ArticleGoogle Scholar
  9. Ranson H, Paton MG, Jensen B, McCarroll L, Vaughan A, Hogan JR, et al. Genetic mapping of genes conferring permethrin resistance in the malaria vector, Anopheles gambiae. Insect Mol Biol. 2004;13(4):379–86.Google Scholar
  10. Feyereisen R. Insect Cytochrome P450. In: Gilbert LI, Latrou K, Gill SS, editors. Comprehensive Molecular Insect Science. 2005Google Scholar
  11. Hudson ME. Sequencing breakthroughs for genomic ecology and evolutionary biology. Mol Ecol Resour. 2008;8(1):3–17.View ArticlePubMedGoogle Scholar
  12. Perry T, Batterham P, Daborn PJ. The biology of insecticidal activity and resistance. Insect Biochem Mol Biol. 2011;41(7):411–22.View ArticlePubMedGoogle Scholar
  13. Georghiou GP. The evolution of resistance to pesticides. Annu Rev Ecol Syst. 1972;3(1):133–68.View ArticleGoogle Scholar
  14. Zhao Q-Y, Wang Y, Kong Y-M, Luo D, Li X, Hao P. Optimizing de novo transcriptome assembly from short-read RNA-Seq data: a comparative study. Bmc Bioinformatics. 2011;12(Suppl 14):S2.View ArticleGoogle Scholar
  15. Haznedaroglu BZ, Reeves D, Rismani-Yazdi H, Peccia J. Optimization of de novo transcriptome assembly from high-throughput short read sequencing data improves functional annotation for non-model organisms. Bmc Bioinformatics. 2012;13:170.PubMed CentralView ArticlePubMedGoogle Scholar
  16. Lin Q, Jin F, Hu Z, Chen H, Yin F, Li Z, et al. Transcriptome analysis of chlorantraniliprole resistance development in the diamondback moth Plutella xylostella. PloS One. 2013;8(8):e72314.PubMed CentralView ArticlePubMedGoogle Scholar
  17. Jackson RJ, Hellen CUT, Pestova TV. The mechanism of eukaryotic translation initiation and principles of its regulation. Nat Rev Mol Cell Biol. 2010;11(2):113–27.PubMed CentralView ArticlePubMedGoogle Scholar
  18. Adolphus, Vanloon PGM, Vaneijk E, Grivell LA. Biosynthesis of the ubiquinol-cytochrome-c reductase complex in yeast - discoordinate synthesis of the 11-kd subunit in response to increased gene copy number. Embo J. 1983;2(10):1765–70.Google Scholar
  19. Keum Y-S. Regulation of Nrf2-Mediated Phase II Detoxification and Anti-oxidant Genes. Biomol Ther. 2012;20(2):144–51.View ArticleGoogle Scholar
  20. Koturbash I, Beland FA, Pogribny IP. Role of microRNAs in the regulation of drug metabolizing and transporting genes and the response to environmental toxicants. Expert Opin Drug Metab Toxicol. 2012;8(5):597–606.View ArticlePubMedGoogle Scholar
  21. Collotta M, Bertazzi PA, Bollati V. Epigenetics and pesticides. Toxicology. 2013;307:35–41.View ArticlePubMedGoogle Scholar
  22. Schnekenburger M, Karius T, Diederich M. Regulation of epigenetic traits of the glutathione S-transferase P1 gene: from detoxification toward cancer prevention and diagnosis. Frontiers in Pharmacology. 2014:5. doi:https://doi.org/10.3389/fphar.2014.00170.
  23. Vandegehuchte MB, Janssen CR. Epigenetics in an ecotoxicological context. Mutat Res-Genet Toxicol Environ Mutagen. 2014;764:36–45.View ArticlePubMedGoogle Scholar
  24. Yu J, Hu S, Ma K, Sun L, Hu H, Zou F, et al. Ribosomal protein S29 regulates metabolic insecticide resistance through binding and degradation of CYP6N3. PloS One. 2014;9(4):e94611.PubMed CentralView ArticlePubMedGoogle Scholar
  25. Oakeshott JG, Farnsworth CA, East PD, Scott C, Han Y, Wu Y, et al. How many genetic options for evolving insecticide resistance in heliothine and spodopteran pests? Pest Manag Sci. 2013;69(8):889–96.PubMed CentralView ArticlePubMedGoogle Scholar
  26. Nascimento ARB, Farias JR, Bernardi D, Horikoshi RJ, Omoto C. Genetic basis of Spodoptera frugiperda (Lepidoptera: Noctuidae) resistance to the chitin synthesis inhibitor lufenuron. Pest Manag Sci. 2015. doi:https://doi.org/10.1002/ps.4057 [Epub ahead of print].
  27. Merzendorfer H, Kim HS, Chaudhari SS, Kumari M, Specht CA, Butcher S, et al. Genomic and proteomic studies on the effects of the insect growth regulator diflubenzuron in the model beetle species Tribolium castaneum. Insect Biochem Mol Biol. 2012;42(4):264–76.View ArticlePubMedGoogle Scholar
  28. Abo-Elghar GE, Fujiyoshi P, Matsumura F. Significance of the sulfonylurea receptor (SUR) as the target of diflubenzuron in chitin synthesis inhibition in Drosophila melanogaster and Blattella germanica. Insect Biochem Mol Biol. 2004;34(8):743–52.View ArticlePubMedGoogle Scholar
  29. Matsumura F. Studies on the action mechanism of benzoylurea insecticides to inhibit the process of chitin synthesis in insects: A review on the status of research activities in the past, the present and the future prospects. Pestic Biochem Physiol. 2010;97(2):133–9.View ArticleGoogle Scholar
  30. Suderman RJ, Dittmer NT, Kanost MR, Kramer KJ. Model reactions for insect cuticle sclerotization: Cross-linking of recombinant cuticular proteins upon their laccase-catalyzed oxidative conjugation with catechols. Insect Biochem Mol Biol. 2006;36(4):353–65.View ArticlePubMedGoogle Scholar
  31. Van Leeuwen T, Demaeght P, Osborne EJ, Dermauw W, Gohlke S, Nauen R, et al. Population bulk segregant mapping uncovers resistance mutations and the mode of action of a chitin synthesis inhibitor in arthropods. Proc Natl Acad Sci U S A. 2012;109(12):4407–12.PubMed CentralView ArticlePubMedGoogle Scholar
  32. Meyer F, Floetenmeyer M, Moussian B. The sulfonylurea receptor Sur is dispensable for chitin synthesis in Drosophila melanogaster embryos. Pest Manag Sci. 2013;69(10):1136–40.View ArticlePubMedGoogle Scholar
  33. Mikolajczyk P, Oberlander H, Silhacek DL, Ishaaya I, Shaaya E. Chitin synthesis in Spodoptera frugiperda wing imaginal disks.1. Chlorfluazuron, diflubenzuron, and teflubenzuron inhibit incorporation but not uptake of C-14 n-acetyl-d-glucosamine. Arch Insect Biochem Physiol. 1994;25(3):245–58.View ArticleGoogle Scholar
  34. Ito Y, Sonobe H. The Role of Ecdysteroid 22-Kinase in the Accumulation of Ecdysteroids in Ovary of Silkworm Bombyx mori. In: Vaudry H, Roubos EW, Coast GM, Vallarino M, editors. Trends in Comparative Endocrinology and Neurobiology vol. 1163. 2009. p. 421–4.Google Scholar
  35. Gangishetti U, Breitenbach S, Zander M, Saheb SK, Mueller U, Schwarz H, et al. Effects of benzoylphenylurea on chitin synthesis and orientation in the cuticle of the Drosophila larva. Eur J Cell Biol. 2009;88(3):167–80.View ArticlePubMedGoogle Scholar
  36. Ikeno T, Ishikawa K, Numata H, Goto SG. Circadian clock gene Clock is involved in the photoperiodic response of the bean bug Riptortus pedestris. Physiol Entomol. 2013;38(2):157–62.View ArticleGoogle Scholar
  37. Huang F-F, Chai C-L, Zhang Z, Liu Z-H, Dai F-Y, Lu C, et al. The UDP-glucosyltransferase multigene family in Bombyx mori. BMC Genomics. 2008;9:563.PubMed CentralView ArticlePubMedGoogle Scholar
  38. Ahn S-J, Vogel H, Heckel DG. Comparative analysis of the UDP-glycosyltransferase multigene family in insects. Insect Biochem Mol Biol. 2012;42(2):133–47.View ArticlePubMedGoogle Scholar
  39. Luque T, O'Reilly DR. Functional and phylogenetic analyses of a putative Drosophila melanogaster UDP-glycosyltransferase gene. Insect Biochem Mol Biol. 2002;32(12):1597–604.View ArticlePubMedGoogle Scholar
  40. Pedra JHF, McIntyre LM, Scharf ME, Pittendrigh BR. Genome-wide transcription profile of field- and laboratory-selected dichlorodiphenyltrichloroethane (DDT)-resistant Drosophila. Proc Natl Acad Sci U S A. 2004;101(18):7034–9.PubMed CentralView ArticlePubMedGoogle Scholar
  41. Silva AX, Jander G, Samaniego H, Ramsey JS, Figueroa CC. Insecticide resistance mechanisms in the green peach aphid Myzus persicae (Hemiptera: Aphididae) I: A transcriptomic survey. PloS One. 2012;7(6):e36366.PubMed CentralView ArticlePubMedGoogle Scholar
  42. Bozzolan F, Siaussat D, Maria A, Durand N, Pottier MA, Chertemps T, et al. Antennal uridine diphosphate (UDP)-glycosyltransferases in a pest insect: diversity and putative function in odorant and xenobiotics clearance. Insect Mol Biol. 2014;23(5):539–49.View ArticlePubMedGoogle Scholar
  43. Brun-Barale A, Hema O, Martin T, Suraporn S, Audant P, Sezutsu H, et al. Multiple P450 genes overexpressed in deltamethrin-resistant strains of Helicoverpa armigera. Pest Manag Sci. 2010;66(8):900–9.PubMedGoogle Scholar
  44. Carvalho RA, Omoto C, Field LM, Williamson MS, Bass C. Investigating the molecular mechanisms of organophosphate and pyrethroid resistance in the fall armyworm Spodoptera frugiperda. PloS One. 2013;8(4):e62268.PubMed CentralView ArticlePubMedGoogle Scholar
  45. Yu SJ, Nguyen SN. Inheritance of carbaryl resistance and microsomal oxidases in the fall armyworm (Lepidoptera: noctuidae). J Econ Entomol. 1994;87(2):301–4.View ArticleGoogle Scholar
  46. Yu SJ. Insecticide resistance in the fall armyworm, Spodoptera frugiperda (J. E. Smith). Pestic Biochem Physiol. 1991;39(1):84–91.View ArticleGoogle Scholar
  47. Pittendrigh B, Aronstein K, Zinkovsky E, Andreev O, Campbell B, Daly J, et al. Cytochrome P450 genes from Helicoverpa armigera: Expression in a pyrethroid-susceptible and -resistant strain. Insect Biochem Mol Biol. 1997;27(6):507–12.View ArticlePubMedGoogle Scholar
  48. Ranasinghe C, Hobbs AA. Isolation and characterization of two cytochrome P450 cDNA clones for CYP6B6 and CYP6B7 from Helicoverpa armigera (Hubner): possible involvement of CYP6B7 in pyrethroid resistance. Insect Biochem Mol Biol. 1998;28(8):571–80.View ArticlePubMedGoogle Scholar
  49. Yang Y, Chen S, Wu S, Yue L, Wu Y. Constitutive overexpression of multiple cytochrome P450 genes associated with pyrethroid resistance in Helicoverpa armigera. J Econ Entomol. 2006;99(5):1784–9.View ArticlePubMedGoogle Scholar
  50. Yang Z, Yang H, He G. Cloning and characterization of two cytochrome P450 CYP6AX1 and CYP6AY1 cDNAs from Nilaparvata lugens Stal (Homoptera : Delphacidae). Arch Insect Biochem Physiol. 2007;64(2):88–99.View ArticlePubMedGoogle Scholar
  51. Zerbino DR, Birney E. Velvet: Algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008;18(5):821–9.PubMed CentralView ArticlePubMedGoogle Scholar
  52. Surget-Groba Y, Montoya-Burgos JI. Optimization of de novo transcriptome assembly from next-generation sequencing data. Genome Res. 2010;20(10):1432–40.PubMed CentralView ArticlePubMedGoogle Scholar
  53. Schulz MH, Zerbino DR, Vingron M, Birney E. Oases: Robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics. 2012;28(8):1086–92.PubMed CentralView ArticlePubMedGoogle Scholar
  54. Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22(13):1658–9.View ArticlePubMedGoogle Scholar
  55. Altschul SF, Madden TL, Schaffer AA, Zhang JH, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.PubMed CentralView ArticlePubMedGoogle Scholar
  56. Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.View ArticlePubMedGoogle Scholar
  57. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36 suppl 1:D480–4.PubMed CentralPubMedGoogle Scholar
  58. Zdobnov EM, Apweiler R. InterProScan - an integration platform for the signature-recognition methods in InterPro. Bioinformatics. 2001;17(9):847–8.View ArticlePubMedGoogle Scholar
  59. Souza TP. Efeito dos inibidores de proteinase de soja no padrão de expressão de proteinases de Spodoptera frugiperda. Piracicaba: Escola Superior de Agricultura Luiz de Queiroz, Universidade de São Paulo; 2013.View ArticleGoogle Scholar
  60. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(T)(−Delta Delta C) method. Methods. 2001;25(4):402–8.View ArticlePubMedGoogle Scholar

Copyright

© Nascimento et al. 2015

Advertisement