Toxicity of Melaleuca alternifolia essential oil to the mitochondrion and NAD⁺/NADH dehydrogenase in Tribolium confusum

EO and chemicals

The EO (density, 0.8978) was purchased from Fujian Senmeida Biological Technology Co., Ltd (Xiamen, China). Terpinen-4-ol (40.09%), γ-terpinene (21.85%), α-terpinene (11.34%), α-terpineol (6.91%), and 1,8-cineole (1.83%) were the major compounds.

Insect culture

A culture of T. confusum was maintained in the laboratory, and the insects were not exposed to any insecticides. For insect culture, the larvae were reared on sterilized whole wheat at 28 ± 1 and 68 ± 5 °C relative humidity under complete darkness. Then, pupae of the same age were collected and transferred to a new container. After emergence, the adults were reared to about 2 weeks of age for use in the subsequent experiments.

Fumigant toxicity assay

The fumigant toxicity of M. alternifolia EO against T. confusum was determined according to our previous protocol (Liao et al., 2016). For oil exposure, 30 adults were exposed to serial dilution doses in sealed gas-tight 300 mL glass jars and incubated for 24, 48, and 72 h at 28 °C. Drops of the oil (1.8, 2.1, 2.5, 3.2, and 4.0 mL) were applied with an Automatic Micro-applicator (Burkard 900X; Burkard Scientific Ltd., Uxbridge, UK) to a piece of filter paper (2 × 3 cm), and the filter paper was attached to the undersurface of the jar lid. Equivalent groups of control adults were treated similarly, but without exposure to the oil. Three biological replicates were maintained for each treatment. For the EO constituents, the protocol for fumigant toxicity was determined using the above-mentioned process, and serial dilutions were prepared and applied to filter paper. In addition, T. confusum specimens exposed to LC₅₀ (6.37 mg/L air) of oil for 12, 24, 36, 48, 60, and 72 h were collected and washed twice or three times with pre-cooled saline, flash-frozen in liquid nitrogen, and stored at −80 °C for the subsequent bioassays.

Transmission electron microscopy of mitochondria

Cells were obtained from the thorax for transmission electron microscopy (TEM) by dissecting the insects at 24, 48, and 72 h after the oil treatment. The samples were fixed in a mixture of 5% glutaraldehyde and 0.1M sodium cacodylate (pH 7.2) for 24 h. After fixation, the samples were washed, dehydrated, and embedded in pure resin, according to the protocol of Correa et al. (2014). After polymerization in gelatin capsules, ultrathin sections were placed on copper grids and subsequently observed and photographed using a transmission electron microscope (HT7700; Hitachi, Tokyo, Japan).

RNA sequencing

Total RNA was extracted from oil treatment and control groups (collected at 24 h) with TRIzol reagent (Kangwei Century Biotechnology Co. Ltd., Beijing, China), according to the manufacturer’s instructions, and treated with DNase I (Sangon Biotech, Shanghai, China). The RNA quality was checked with a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Library construction and Illumina sequencing were performed at BGI-Tech (Wuhan, China). For cDNA library construction, five μg of RNA per sample from three biological replicates were combined and used. Two cDNA libraries were constructed for the oil treatment and control groups. For Illumina sequencing, which followed the protocol of the Illumina TruSeq RNA Sample Preparation Kit (BGI-Tech, Wuhan, China), 2 × 100 bp paired-end reads were sequenced using Illumina HiSeq™ 4000 (Illumina Inc., San Diego, CA, USA), with the depth of six g for each sample. The reads were submitted to the NCBI Sequence Read Archive (SRA; accession number, SRS2593554).

Bioinformatic analyses

The reads for the treatment and control groups were mapped to the 165.944-Mb T. castaneum reference genome obtained from NCBI (BioProjects: PRJNA12540) by using TopHat v.2.08 (Kim et al., 2013a), with quality aware alignment algorithms (Bowtie v.2.2.5) (Langmead et al., 2009).

The raw RNA-seq reads were assessed for quality with FastQC (version 0.11.4; Babraham Bioinformatics, Cambridge, UK) and saved as FASTQ files with default parameters (Cock et al., 2010). Then, de novo assembly of the clean reads was performed using the Trinity method (version 2.0.6) (Grabherr et al., 2011). All the unique Trinity contigs were analyzed using BlastX (E-value < 10⁻⁵) against the protein databases Nr (Agarwala et al., 2016), Nt (Agarwala et al., 2016), COG (Tatusov et al., 2000), KEGG (Kanehisa & Goto, 2000), Swiss-Prot, and InterPro using InterProScan5 with default parameters. To annotate the assembled sequences with GO terms, Nr Blast results were imported into Blast2GO (Conesa et al., 2005).

Transcript abundance was calculated as fragments per kilobase of transcript per million fragments mapped (FPKM) for each sample (Li & Dewey, 2011). Differential gene expression analysis (fold changes) and related statistical significance in pair-wise comparison were performed using the DESeq program (http://www-huber.embl.de/users/anders/DESeq/) (Anders & Huber, 2010). The differentially expressed genes (DEGs) were identified using a false discovery rate (FDR) threshold ≤ 0.001 and absolute value of log2Ratio ≥ 1 (Hao et al., 2016). Genes with an adjusted p-value were used for controlling FDR, and those with a threshold < 0.05 were classified as differentially expressed (Ma et al., 2015).

For each DEG, GO and KEGG enrichment analyses were conducted using the DESeq R package (http://www.geneontology.org/ and http://www.genome.jp/kegg/, respectively). The GOslim annotations results were then classified into three main classes: molecular function, biological process, and cellular component. The KEGG database was used to identify significantly enriched metabolic pathways or signal transduction pathways.

Quantitative real-time PCR

Quantitative real-time (qRT-PCR) was used to further validate and quantify the RNA levels for 20 selected genes that encode NADH or NAD⁺ by using the iCycler iQ Real-time Detection System (Bio-Rad, Hercules, CA, USA). Gene-specific primers were designed using Primer Premier 5, and the sequences are listed in Table S1. The house-keeping gene glyceraldehyde 3-phosphate dehydrogenase was used as the reference gene, as proposed by Pan et al. (2015). For the qRT-PCR analysis, cDNA templates were diluted 20-fold in nuclease-free water. Then, mRNA levels were measured in triplicate (technical repeats) with qPCR by using the SYBR Green Master Mix (Vazyme Biotech Co., Ltd, Nanjing, China), according to the manufacturer’s instructions. PCR amplification was performed in a total volume of 20.0 μL containing 10.0 μL of the SYBR Master Mix, 0.4 μL of each primer (10 μM), 2.0 μL of cDNA, and 7.2 μL of RNase-free water. The amplification procedure was composed of an initial denaturation step at 95 °C for 5 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s and the melting curve step at 95 °C for 15 s, 60 °C for 60 s, and 95 °C for 15 s. Gene expression was quantified (mean ± SD) as relative fold change by using the 2^−ΔΔCT method (Schmittgen & Livak, 2008).

Measurement of intracellular NAD⁺/NADH ratio

Both oxidized and reduced forms of intracellular NAD were determined using an NAD(H) quantification kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). Briefly, 0.1 g of the test insects were collected at 12, 24, 36, 48, 60, and 72 h and extracted with one mL of NAD⁺/NADH extraction buffer in three freeze/thaw cycles. The samples were centrifuged at 10,000×g for 5 min at 4 °C. Then, 0.5 mL of the extracted NADH or NAD⁺ supernatant was transferred to a centrifuge tube and neutralized with an equal volume of the opposite extraction buffer. The samples were centrifuged at 10,000×g for 10 min at 4 °C and then used for the subsequent bioassays. NADH or NAD⁺ cycling mix was prepared according to the manufacturer’s protocol (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). Finally, absorbance was measured at 570 nm. In addition, the concentration of the total protein was determined using the total protein quantitative assay (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). Three replicates were used for each treatment, and each replicate was determined three times.

Statistical analysis

The mortality rates observed in the toxicity bioassays were corrected for the control group by using Abbott’s formula (Abbott, 1925). All data are expressed as mean ± SE values of three independent experiments and analyzed using one-way nested analysis of variance and unpaired sample t-test. A significant difference was accepted at a p-value < 0.05. An extremely significant difference was accepted at p-value < 0.01. The LC₅₀ values were evaluated using probit analysis (Hubert & Carter, 1990), and corresponding confidence intervals at 95% probability were obtained using IBM SPSS Statistics 22.0 (IBM, Armonk, NY, USA). Figures depicting the effects of the EO on enzymatic activities and the qRT-PCR results were created using Origin Pro 9.0 (Origin Lab Corporation, Northampton, MA, USA).

Fumigant toxicity of M. alternifolia EO and constituents

To investigate the toxicity of the M. alternifolia EO against T. confusum adults, we performed the fumigation assay. The results show that M. alternifolia EO has potent fumigant toxicity (Fig. 1A), and the effect of fumigation gradually increased over time (24, 48, and 72 h); the corresponding LC₅₀ values were 7.45, 7.09, and 6.37 mg/L air, respectively (Fig. 1B). The largest dose of 11.97 mg/L air EO caused 91.11%, 97.78%, and 98.86% mortality, respectively, in the T. confusum adults.

In particular, terpinen-4-ol was the most potent toxicant with an LC₅₀ value of 3.83 mg/L air (Fig. 1C). In the M. alternifolia EO, terpinen-4-ol was the main component (40.09% of the EO), indicating that terpinen-4-ol is the major contributor to the fumigant toxicity of the EO. In addition, γ-terpinene and α-terpinene exhibited weaker fumigant toxicity (LC₅₀ = 28.52 and 44.53 mg/L air, respectively) against T. confusum.

TEM of mitochondria

An ultra-structural examination of the morphology of the mitochondria from untreated and oil-fumigated T. confusum larvae is shown in Fig. 2. In the untreated T. confusum larvae, the mitochondria have highly electron-dense cristae, membranes, and matrix (Figs. 2A, 2C and 2E). However, the mitochondria in the columnar and regenerative nidi cells from the thorax of the oil-treated T. confusum larvae had undergone ultra-structural changes detected by the vacuolization of the mitochondrial matrix (Figs. 2B, 2D and 2F), when compared with the non-fumigated adults. The vacuolization increased with time after the oil treatment and, in severe cases, caused fragmentation of the mitochondria.

Illumina sequencing and de novo assembly

To obtain a global, comprehensive overview of the T. confusum transcriptome, RNA was extracted from the treatments and control groups. A total of 126,280,032 paired-end reads (100 bp) were generated from the samples by using the Illumina HiSeq™ 4000 platform. Then, 89,342,546 clean reads were obtained by preprocessing and filtering the reads (low-quality sequences were removed; Table 1). Subsequently, the clean reads were subjected to transcriptome assembly by using the Trinity software package (Grabherr et al., 2011), and 28,885 assembled unigenes were generated using overlapping information from high-quality reads, which accounted for 36,998,010 bp (Table 1). Of the assembled unigenes, approximately 38.54% were ≤600 bp and 61.46% were >500 bp. The average length of the unigenes was 1,280 bp, with an N50 length of 2,097 bp and mean length of 1,280 bp. The length distribution of the unigenes is shown in Fig. 3A.

Functional annotation of T. confusum transcripts

All the assembled unigenes were aligned against seven public databases (Table 1). Of the 28,885 assembled unigenes, 23,160 (80.18%) exhibited sequence similarity to a sequence within the Nr database; 23,571 (81.60%) unigenes were annotated in at least one database, indicating that just a few unigenes (18.40%) could not be identified. The homologous genes that showed the best match (54.24%) were from T. castaneum (91.82%). On the basis of the Nr annotation, GO functional analysis of the unigenes was performed. A total of 6,333 (21.92%) unigenes were assigned to the biological process, molecular function, and cellular component categories, including 57 GO terms (Fig. 3B). In addition, 18,074 (62.57%) unigenes were divided into 42 subcategories and 295 KEGG pathways by using the KEGG annotation system with default parameters to predict the metabolic pathways (Fig. 3C).

Differential expression analysis and pathway enrichment

The sequence analysis and annotation of all the unigenes in S. zeamais fumigated by the M. alternifolia EO provided some valuable information for analyzing the T. confusum transcriptome. From the 23,571 unigenes identified in the analysis, we chose to focus on transcripts encoded by the genes associated with known mechanisms to cope with xenobiotic compounds, including quantitative or qualitative changes in major detoxification enzymes and transporters to decrease exposure (pharmacokinetic mechanisms) or changes in target site sensitivity (pharmacodynamic mechanisms) (Bajda et al., 2015). Specifically, changes in the expression levels of four classes of enzymes and proteins (GST, CarE, cytochrome P450 monooxygenases, and mitochondrial respiratory chain-related proteins) were investigated to determine whether patterns emerged in the upregulation or downregulation of specific transcripts. The transcriptome of T. confusum showed that the largest and most abundant group was ATPase transporters, followed by cytochrome P450s; some of them may be involved in insecticidal mechanisms. The transcriptome also showed five possible NAD⁺/NADH dehydrogenase transcripts, which may be the main targets for the EO.

For comparison, FPKM of each transcript was calculated to estimate the expression levels between the oil-fumigated and oil-free samples. The important DEGs (999 upregulated and 1,209 downregulated) were identified on the basis of threshold FDR < 0.01 and fold change 2 between the oil-fumigated and oil-free samples. To annotate these DEGs, both GO and KEGG functional analyses were performed.

The GO annotation analysis classified 632 DEGs into three GO categories and 339 terms (Fig. S1A). In the molecular function category, 560 DEGs were classified into 11 terms, namely, antioxidant activity, binding, catalytic activity, electron carrier activity, enzyme regulator activity, guanyl-nucleotide exchange factor activity, molecular transducer activity, nucleic acid binding transcription factor activity, receptor activity, structural molecule activity, and transporter activity.

Among the DEGs, 1,180 unigenes were mapped to 287 different KEGG pathways and five categories (Fig. S1B). According to the threshold of Q-value < 0.05, 22 pathways were significantly enriched (Table S2). Many DEGs were significantly enriched in the metabolism pathways associated with respiration and metabolism of xenobiotics, suggesting that abnormal respiration and metabolic disorders occurred in the T. confusum adults after fumigation with the M. alternifolia EO. In addition, 92 possible insect hormone biosynthesis transcripts, some of which are known targets of chlorbenzuron, were detected (Xu, Xu & Wu, 2017).

To verify the expression patterns of the DEGs involved in metabolism, 20 genes were selected for qRT-PCR analysis. As shown in Fig. 4, similar trends of upregulation/downregulation of the selected DEGs were observed between the qRT-PCR and transcriptome data, confirming the accuracy of our transcriptome profiling.

NAD⁺/NADH ratio in T. confusum fumigated with the M. alternifolia EO

On the basis of a previous study on S. zeamais and the above-mentioned results, the NAD⁺/NADH ratio in T. confusum fumigated with the M. alternifolia EO was measured to investigate whether the EO acts on NAD⁺/NADH. In the non-fumigated insects, a decrease in NAD⁺ and NADH levels was observed over the course of 24–48/60 h, which may be affected by starvation. Further, we found that treatment with 6.37 mg/L EO significantly increased NAD⁺ (Fig. 5A) but decreased NADH (Fig. 5B) levels at 12–48 h, when compared with the non-fumigated samples; however, the opposite trend was observed after 60 h. The ratio of NAD⁺/NADH in T. confusum from 12 to 60 h after treatment decreased (significantly in 24–48 h) and increased after 60 h, but not effectively (Fig. 5C).