Gut transcriptome of two bark beetle species stimulated with the same kairomones reveals molecular differences in detoxification pathways.

Dendroctonus bark beetles are the most destructive agents in coniferous forests. These beetles come into contact with the toxic compounds of their host's chemical defenses throughout their life cycle, some of which are also used by the insects as kairomones to select their host trees during the colonization process. However, little is known about the molecular mechanisms by which the insects counteract the toxicity of these compounds. Here, two sibling species of bark beetles, D. valens and D. rhizophagus, were stimulated with vapors of a blend of their main kairomones (α-pinene, β-pinene and 3-carene), in order to compare the transcriptional response of their gut. A total of 48 180 unigenes were identified in D. valens and 43 704 in D. rhizophagus, in response to kairomones blend. The analysis of differential gene expression showed a transcriptional response in D. valens (739 unigenes, 0.58-10.36 Log2FC) related to digestive process and in D. rhizophagus (322 unigenes 0.87-13.08 Log2FC) related to xenobiotics metabolism. The expression profiles of detoxification genes mainly evidenced the up-regulation of COEs and GSTs in D. valens, and the up-regulation of P450s in D. rhizophagus. Results suggest that terpenes metabolism comes accompanied by an integral hormetic response, result of compensatory mechanisms, including the activation of other metabolic pathways, to ensure the supply of energy and the survival of organisms which is specific for each species, according to its life history and ecological strategy.

Introduction

Dendroctonus bark beetles (Coleoptera: Scolytidae) are the most destructive agents of coniferous forests because their populations produce large-scale outbreaks causing the death of hundreds of thousands of trees, resulting in ecological disturbances, changes in the landscape, as well as significant economic losses [1]. These bark beetles are conservative in the use of hosts, because they colonize only trees of the genera Larix, Pseudotsugae, Picea and Pinus (Pinaceae), mainly pines and spruces, throughout their distribution range [2].

For the selection of their host trees, pioneer beetles use volatile terpenes released by the trees as primary attractants (kairomones). These attractants, along with the releasing of de novo synthesized pheromones and compounds derived from terpenes metabolism by these beetles, favor massive attacks on trees and their successful colonization [3], [4]. Once in the tree, they must overcome the defensive mechanisms of host trees, mainly chemical defenses integrated by monoterpenes, non-volatile diterpenes, sesquiterpenes, and phenolic compounds present both in the constitutive and induced resin [5], [6]. These compounds come in direct contact with beetles from the moment they drill into the tree’s bark, and their toxicity produces severe damage to organelles, cells, or morpho-physiological systems [7], [8], affecting pheromone production and even causing the death of insects and their symbiotic microorganisms [9], [10]. Unfortunately, little is known about molecular mechanisms to counteract the toxicity of these compounds.

The metabolism of xenobiotics activates the detoxification process in phytophagous insects. Detoxification is a complex process of chemical transformations involving several multigenic families performing different types of reactions and excretion processes: 1) functionalization (phase I), 2) conjugation (phase II), and transportation (phase III) [11]. Enzymes associated with functionalization, such as cytochrome P450 monooxygenases (CYPs or P450s) and carboxylesterases (COEs), catalyze reactions that introduce reactive and polar groups into their substrates via oxidation, reduction, or hydrolysis reactions to facilitate sequestration or excretion into more water-soluble compounds. Conjugation enzymes such as glutathione S-transferases (GSTs) conjugate or convert lipophilic compounds into more hydrophilic compounds [12]. Both types of reactions and enzymes occur in all different developmental stages, sex, and time, and taken place simultaneously or independently in different tissues, organs and morpho-physiological systems [13], [14]. Transportation processes are carry out mainly by ABC transporters enzymes that hydrolyze ATP and move hydrophobic or conjugated products out of cells for excretion [15].

The Dendroctonus spp. gut is a system that plays a vital role in digestive, detoxification and pheromone synthesis processes [16], [17], [3]. In particular, the integration of detoxification enzymes into specific physiological pathways in the gut, related to pheromone production, might be a crucial part of the metabolic exaptation from these bark beetles throughout their evolutionary history. These enzymes are related to the beetles’ ability to metabolize, sequester and/or solubilize tree’s resin compounds within a short time in a changing environment, but also some of them participate simultaneously or later in the pheromone synthesis [4].

Several studies in bark beetles have characterized and analyzed changes in the expression level of specific genes encoding enzymes such as P450s, COEs, and GSTs under different experimental conditions and stimuli with terpenes (e.g. insects stimulated with terpene vapors in the lab., phloem-fed insects in the lab and field, expression in different tissues, developmental stages, exposure times, sexes, and colonization stages) using different techniques [18], [19], [20], [21], [22]. These studies showed direct upward and downward relationships in gene expression and the applied stimulus under the different experimental conditions assayed, as well as their direct participation, mainly of specific P450s, in the metabolism of terpenes [23], [24], [25].

To get an overview of genes involved in different metabolic processes, integral and specific tissue transcriptoms have been recognized as a powerful approach [26]. These studies have been conducted in single species of Dendroctonus, which have demonstrated the molecular richness associated with chemoreception, pheromone biosynthesis and terpenes detoxification [27], [28], [29], [30], [31], [32]. However, comparative studies between closely related species with different biological, behavioral, and ecological attributes, but with similar chemical communication systems, have not been performed. Dendroctonus valens LeConte and D. rhizophagus Thomas and Bright are sibling species that differentiated during the late Pleistocene [33], [34]. Both species are attracted by the same kairomones of pine trees: α-pinene, β-pinene and 3-carene, and produce the same oxygenated-monoterpenes: trans-verbenol, cis-verbenol, myrtenal, mirtenol, and verbenone [35], [36], [37]. These suggest that they have similar chemical communication systems, excepting because D. valens produces de novo the pheromones frontalin and exo-brevicomin [38], [39].

Several studies have documented the significant transcriptional activity of specific genes in the gut of D. valens and D. rhizophagus under different experimental conditions [40], [41], [19], [21]. In this research, a comprehensive gut-specific transcriptome was generated, with the aim of knowing the global transcriptional activity occurring in this system after exposing insects of both species to vapors of a blend of their main kairomones: α-pinene, β-pinene and 3-carene. We also documented the diversity of genes P450s, COEs and GSTs associated with the detoxification process (phase I and II) and compared the intra- and interspecific expression profiles of these genes.

Materials and methods

Insect collection and treatment

Pre-emerged adults of D. valens and D. rhizophagus were collected directly from naturally infested Arizona pines (P. arizonica Engelm.) in the Sierra del Tigre, Zapotlán El Grande Municipality, Jalisco State (19°42′10″ N, 103°27′45″ W; 1520 m a.s.l.), and in San Juanito, Bocoyna Municipality, Chihuahua State (29°44′17″ N, 108°15′43.7″ W; 2288 m a.s.l.) Mexico, respectively. The sex of beetles was determined based on the shape of the seventh abdominal tergite and stridulation by males.

Since beetles naturally come in contact with these compounds via the cuticle, as well as the respiratory and the digestive system, the insects were stimulated with vapors of a blend (1:1:1) from α-pinene, β-pinene and 3-carene (Synergy Semiochemical, Canada) to avoid as much as possible the stimuli by contact or feeding. A group of ten insects per species, with a sexual ratio of 1:1, were placed inside Petri dishes (100 mm × 15 mm). A piece of filter paper impregnated with 100 µL (0.5 mM) of the kairomones blend was adhered to the inner surface of the lid of the Petri dishes, preventing the insects from getting into physical contact with the chemical compounds. To minimize the loss of compounds by evaporation, Petri dishes were closed and sealed with Parafilm (Bemis Company Inc.). Insects were stimulated during 8 h in darkness. Ten non-stimulated insects of each bark beetle constituted the control groups, and two replicates were realized for each species.

Beetles were killed by immersion in Phosphate-Buffered-Saline (PBS 137 mmol/L NaCl, 2.7 mmol/L KCl, 10 mmol/L Na2HPO4, 2 mmol/L KH2PO4 at pH 7.4). Using a stereoscopic microscope and fine-tipped tweezers, the gut of the insects was extracted. Guts corresponding to each replicate per species and control group were pooled separately into Eppendorf tubes containing 0.2 mL of TRI Reagent® solution (Ambion® by Life Technologies, Carlsbad, CA, USA) and homogenized with a sterile pestle. All tubes were filled with TRI Reagent® solution up to 1 mL and stored at −80 °C until total RNA extraction.

Total RNA isolation, library construction, and Illumina sequencing

The total RNA of each sample was extracted using RiboPureTM RNA Purification kit (Ambion® by Life Technologies, Carlsbad, CA, USA), following the manufacturer’s protocol. RNA quality and quantity were measured using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA). From each sample of total RNA, cDNA was synthesized using TruSeq RNA Sample Preparation Kit v2 (Illumina Inc., San Diego, CA, USA) according to the manufacturer’s instructions. The cDNA was later sequenced using Illumina HiSeq 2000 platform (Illumina Inc., San Diego, CA, USA), with a maximum read length of 2 × 100 bp. The construction of libraries and their sequencing were performed at the National Laboratory of Genomics for Biodiversity, LANGEBIO, CINVESTAV-Irapuato, Mexico.

Transcriptome assembly and functional annotation

The quality of the raw readings from each library was checked and visualized with FastQC v.0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The Illumina adapters and low-quality reads were trimmed with Trimmomatic v.0.36 [42], using a sliding window of kmers = 5 to remove leading or trailing bases with an average of Phred quality score lower than 28, as well as reads below 25 pb. The cleaned reads corresponding to each species were used to perform the de novo transcriptome assemblies with the short-read assembly program Trinity v.2.8.5 [43], using default parameters. Transcript redundancies were removed by clustering with an identity threshold of 95% in CD-HIT v.4.6.6 [44], the longest transcripts were considered as representative for each cluster.

The multi-level quality evaluation of both transcriptomes was achieved in three steps: (1) Evaluation of the percentage of properly paired reads mapped to the final assemblies with Bowtie2 v.2.2.9. [45] using the parameters suggested by Trinity; (2) Determination of the percentage of orthologous genes from endopterygota_odb10 database present in each assembly, using BUSCO v.5 [46]; and (3) Estimation of the number of full-length transcripts against complete genes from the genome of Dendroctonus ponderosae by the Perl script “analyze_blastPlus_topHit_coverage.pl” of Trinity.

Using the Trinotate v.3.0.0 pipeline (https://trinotate.github.io/), transcriptomes from D. valens and D. rhizophagus were annotated. BLASTX with E-value cut-off ≤ 10−5 was used for homology search of all putative transcripts with the Swissprot-Uniprot database (https://www.Uniprot.org). Open reading frames (ORFs) of putative transcripts were predicted using TransDecoder v.3.0.1 (https://transdecoder.github.io/), and only predicted ORFs that were at least 100 amino acids long, whether partial or complete, were selected and subsequently identified using BLASTP (E-value cut-off ≤ 10−5). Identification of conserved protein domains was achieved using Hammer v.3.1b1 [47] in the Pfam domain database (pfam.xfam.org). To determine the distribution of putative transcripts functions of each assembly, the results of the Gene Ontology (GO) were loaded into WEGO v.2.0 [48]. The WebMGA [49] was employed to identify orthologous and paralogous eukaryotic proteins in the KOG database (EuKaryotic Orthologous Groups). Lastly, the assignment of putative transcripts into metabolic pathways of the KEGG database (Kyoto Encyclopedia of Genes and Genomes) was performed on the GhostKOALA server (https://www.kegg.jp/ghostkoala).

Differential expression analysis

Through the use of Bowtie2 software and allowing not more than one nucleotide mismatch, clean reads of each cDNA library were mapped on their corresponding transcriptome. Then, RSEM v.1.3.2 [50] was used to identify the number of unambiguous reads, and to estimate the expression levels of transcripts. The parameters recommended by Trinity were used in Bowtie2 and RSEM. The TMM (Trimmed Mean of M values) method was used to normalize the abundance and lowly expressed unigenes < 5 TPM (transcripts per million reads) were filtered. Differentially expressed unigenes (DEGs) in the samples of kairomone-stimulated insects versus control of each species were identified using the edgeR with a log2 fold change (Log2FC) > 0.5 (up-regulated) or Log2FC < − 0.5 (down-regulated) and false discovery rate (FDR) < 0.05.

GO and KEGG enrichment analysis of DEGs

In order to predict the differentially expressed functions between samples with kairomones or control, an analysis of enrichment of GO terms and KEGG pathways was realized with clusterProfiler v.3.10.1 [51]. A hypergeometric test is used by this package to correct the p-value and control the false discovery rate for each GO term/ KEGG pathway. Significant GO terms and KEGG pathways were selected with p-value < 0.05 and q-value < 0.05.

Gene detoxification diversity

In the annotation report generated with Trinotate for each species, genes related to the detoxification process were searched based on full names, abbreviated names, and protein domains of the enzymatic families P450s (PF00067), COEs (PF00135), and GSTs (PF00043, 02798 and 13417). The identity of putative genes, greater than or equal to 250 codons for P450s and COEs, and greater than or equal to 100 codons for GSTs, was verified in BLASTP (E-value cutoff <10−5) against the NCBI-nr database.

Characterization in silico

A protein isoelectric point calculator server (https://isoelectric.org) was used to predict physicochemical characteristics such as molecular weight and the isoelectric point (pI) of complete detoxification of enzymes. Subcellular localization was predicted with Protcomp-AN v.9.0 (Softberry, Inc.). Additionally, conserved motifs of each family were detected by Multiple Em for Motif Elicitation (MEME) suite v.5.3.3 (https://meme-suite.org/meme/tools/meme) with the following parameters: mod = zoops, nmotifs = 30, minw = 6, and maxw = 40.

Phylogenetic analysis

To analyze the phylogenetic relationships between the putative detoxification genes found in D. valens and D. rhizophagus, a search for P450s, COEs and GSTs amino acid sequences annotated into the genomes from Dendroctonus ponderosae (GCA_000355655.1), Anoplophora glabripennis (GCF_000390285.2,), and Leptinotarsa decemlineata (GCF_000500325.1) [52], [53], [54] deposited in the GenBank (https://www.ncbi.nlm.nih.gov) was performed. Additionally, the sequences obtained were concatenated in a single file with which a BLASTP analysis was performed, with an E-value cut-off ≤ 10−5 for genomes used, as well as for the transcriptome from D. frontalis (GAFI00000000) [55], in order to recover all the present sequences. The sequences obtained from the transcriptomes of D. valens and D. rhizophagus, as well as from the genomes of other beetles and the sequences from D. armandi deposited in the GenBank (AHL26982, AHL26983, ALD15890-ALD15933 P450s; AYN64423.1- AYN64430.1 COEs; KJ783317, KJ783316, KJ637332, KP258217-KP258224 GSTs) were used. Multiple alignments were carried out using MUSCLE v.3.8.31 [56], and the phylogenetic tree was performed with PhyML v.3.0 in ATGC Montpellier Bioinformatic Platform (https://www.atgc-montpellier.fr/phyml). The consistency at each node was assessed by aLRT SH-like and bootstrapping after 100 pseudoreplicates. The best evolutionary model for protein data was selected using SMS software and follow the Akaike Information Criterion in the same platform.

Results

About 40 million high-quality paired-end reads were combined to assemble de novo the transcriptome for each bark beetle species. A total of 48 180 unigenes were obtained in D. valens and 43 704 in D. rhizophagus, with a N50 of 2462 and 2701 pb, respectively (Table 1). The size analysis showed that 21.08% and 23.05% of unigenes from D. valens and D. rhizophagus, respectively, are in the range from 1 001 to 3 000 nt (Supplementary Fig. 1A).

Table 1

Statistical summary of the de novo gut transcriptome analysis from D. valens and D. rhizophagus.

	D. valens				D. rhizophagus
	Control		With kairomones		Control		With kairomones
	Dval-C1	Dval-C2	Dval-K1	Dval-K2	Drhi-C1	Drhi-C2	Drhi-K1	Drhi-K2
Raw reads	12 592 430	13 409 535	12 476 737	10 471 278	12 203 329	13 216 638	12 961 902	12 740 627
Clean reads	10 607 553	11 467 453	10 581 596	8 767 541	10 031 115	11 238 692	10 998 277	10 628 661
Q30%	86.91	88.34	86.91	86.91	86.91	86.91	86.91	88.34
GC%	42	43	43	43	43	43	43	44
Transcripts	59 538				55 993
N50	2829				3042
Median conting length (bp)	509				627
Unigenes	48 180				43 704
Cluster of unigenes	7730				8153
Singleton of unigenes	40 450				35 551
N50	2462				2701
Median conting length (pb)	427				487
Annotated in Swiss-Prot	17 814 (36.97%)				15 995 (36.60%)
Annotated in Pfam	14 398 (29.88%)				7711 (17.64%)
Annotated in GO	17 268 (35.84%)				15 611 (35.72%)
Annotated in KOG	12 964 (26.91%)				12 845 (29.39%)
Annotated in KEGG	11 160 (23.16%)				11 057 (25.30%)
Annotated in all database	9453 (19.62%)				5708 (13.06%)
Annotated in at least one database	19 166 (39.78%)				16 695 (38.20%)

Completeness analysis of de novo assemblies showed that >97% clean reads were properly matched in the assembly of each species. The BUSCO evaluation showed that of the 2124 genes included in the endopterygota_odb10 database, the 92.01% and 92.18% were found in the transcriptomes from D. valens and D. rhizophagus as complete genes; only the 3.07% and 2.74% fragmented, and 4.92% and 5.08% missing (Supplementary Fig. 1B). The unigenes identified in D. valens (75.28%) and D. rhizophagus (74.97%) had the best alignment hit at 70% coverage, with proteins identified in the genome of D. ponderosae (E-value of 1e−20).

From the 48 180 unigenes recorded in D. valens and the 43 704 in D. rhizophagus, 19 166 (39.79%) and 16 695 (38.20%), respectively, were annotated in at least one database (Swiss-Prot, Pfam, GO, KOG, and KEGG; Table 1). The unigenes from D. valens and D. rhizophagus were mainly assigned to the following functional categories and dominant subcategories in GO terms: Cellular Components at “cell” and “cell part”, Molecular Function at “binding” and “catalytic activity”, and Biological Process at “cellular process” and “metabolic process” (Fig. 1A).

Fig. 1

D. valens (dark colors) and D. rhizophagus (light colors) unigenes functional annotation. A) GO clasification. The GO terms were classified into of cellular component, molecular function, and biological process. B) KOG classification in four Macrogroups and 25 categories. Cellular processes and signaling: D, cell cycle control cell division, chromosome partitioning; M, cell wall/membrane/envelope biogenesis; N, cell motility; O, posttranslational modification, protein turnover, chaperones; T, signal transduction mechanisms; U, intracellular trafficking, secretion, and vesicular transport; V, defense mechanisms; W, extracelular structures; Y, nuclear structure; Z, cytoskeleton. Poorly characterized: R, general function prediction only; S, function unknown. Infotmation Storage and Processing: A, RNA processing and modification; B, chromatin structure and dynamics; J, translation, ribosomal structure and biogenesis; K, transcription; L, replication, recombination and repair. Metabolisms: C, energy production and conversion; E, amino acid transport and metabolism; F, nucleotide transport and metabolism; G, carbohydrate transport and metabolism; H, coenzyme transport and metabolism; I, lipid transport and metabolism; P, inorganic ion transport and metabolism; Q, secondary metabolites biosynthesis, transport and catabolism. C) KEGG classification in five major metabolic categories: M, metabolism; G, genetic information processing; E, environmental information processing; C, cellular processes; O, organismal systems.

The categorization based on the KOG database assigned the unigenes of both species to 25 groups (A-Z) into four main macro-groups: Cellular Process and Signaling at “signal transduction mechanisms (T)”, Information Storage and Processing at “transcription (K)”, Metabolism at “lipids transport and metabolism (I)”, and Poorly Characterized at “general function prediction only (R)” (Fig. 1B).

The KEGG database assigned the unigenes of both bark beetles into 396 and 392 pathways. The best represented metabolic pathways were: Metabolism at “global and overview maps”; Genetic Information Processing at “translation”, Environmental Information Processing at “signal transduction”, Cellular Processes at “cell growth and death”, and Organismal Systems at “endocrine system” (Fig. 1C).

Overall expression pattern

Based on a threshold of 0.5 Log2FC and FDR < 0.05, only the 2.39% (1 152) of the total unigenes were differentially regulated in stimulated insects of D. valens versus control, whereas only 1.35% (592 unigenes) were differentially expressed in D. rhizophagus (Table 2, Supplementary Table S1).

Table 2

Number of differentially regulated unigenes.

	Total	Control a	With kairomones a
D. valens	1152	413	739
D. rhizophagus	592	270	322

Log2FC > 0.5 and FDR < 0.05.

The top 20 of DEGs in the control of both species correspond to predicted enzymes related to basic cell functions, while the top 20 in insect stimulated with blend of kairomones corresponded to predicted enzymes related to digestive and detoxification functions. Interestingly, the specific functions were different between bark beetles (Supplementary Table S2).

The enrichment analysis of up-regulated unigenes in response to kairomones (Supplementary Table S3) showed that the most significant GO terms for D. valens were associated with: “Microtubule” (GO:0005874, 29 unigenes), “Carbohydrates metabolic process” (GO:0005975, 24 unigenes), “Lysosome” (GO:0005764, 21 unigenes), and “Digestion” (GO:0007586, 17 unigenes) (Fig. 2A). Only seven significant KEGG pathways were found, being the most important: “Glutathione metabolism” (K11140, four unigenes) and “ABC transporters” (K05673, five unigenes) (Fig. 2C). For D. rhizophagus, the most significant GO terms were: “Cell surface” (GO:0009986, 13 unigenes), “Heme binding” (GO:0020037, 12 unigenes), “Protein homotetramerization” (GO:0051289, seven unigenes), and “Ion transport” (GO:0006811, seven unigenes) (Fig. 2B). Eighteen KEGG pathways were significant in this species, of which the most important were: “Choline metabolism in cancer” (K08202, six unigenes); “Cytochrome P450” (K14999, five unigenes); and “Oxidative phosphorylation” (K11352, four unigenes) (Fig. 2D).

Fig. 2

Main GO terms and KEGG pathway enriched in response to the blend kairomone. (a) Ration of GO terms in D. valens. (b) Ration of GO terms in D. rhizophagus. (c) Ration of KEGG pathways in D. valens. (d) Ration of KEGG pathways in D. rhizophagus. The color scale represent enrichment significance p < 0.05.

Detoxification genes

Our results show a minimal difference in the total number of detoxification genes expressed between D. valens (117) and D. rhizophagus (118) (Table 3). Yet, only some of them were up-regulated in these bark beetles (Table 4). These genes represent around 73% of the diversity of P450s, COEs, and GTSs present in genomes from D. ponderosae (159), 57% L. decemlineata (202), and 43% A. glabiprennis (272) (Table 3, Supplementary Table S4).

Table 3

Comparison of detoxification genes in different coleoptera.

	Curculionidae					Cerambicidae	Chrysomelidae
	Dendroctonusvalensa	Dendroctonusrhizophagusa	Dendroctonusfrontalisa	Dendroctonusarmandib	Dendroctonusponderosaec	Anoplophraglabripennisc	Leptinotarsadecemlineatac
P450 clans
CYP mitochondrial	6	5	2	6	8	13	14
CYP2	5	6	0	5	7	7	9
CYP3	27	25	20	19	40	71	47
CYP4	14	15	6	15	21	37	20
Total	52	51	27	46	76	128	90
COE clades
Dietary and detoxification functions
A clade	21	22	18	6	24	46	41
C clade	6	6	6	0	11	8	15
Hormone and semiochemical processing
D clade	1	1	1	0	1	6	8
E clade	3	2	3	1	3	33	4
F clade	1	1	0	1	1	2	1
Neurodevelopmental functions
H clade	1	1	0	0	1	1	1
I clade	0	0	0	0	1	1	1
J clade	1	1	0	0	2	2	2
K clade	1	1	1	0	1	1	1
L clade	1	1	1	0	5	6	5
M clade	1	1	0	0	1	1	1
Total	37	37	30	8	51	107	80
Cytosolic GST classes
Zeta	1	1	1	0	1	1	1
Theta	1	1	0	1	2	2	2
Omega	3	3	2	2	3	4	6
Sigma	6	7	6	2	6	12	3
Delta	2	3	2	0	3	7	3
Epsilon	12	12	9	6	15	10	13
Others	3	3	3	0	2	1	4
Total	28	30	23	11	32	37	32

Detoxification genes identified in the transcriptomes of D. valens (this study), D. rhizophagus (this study), and D. frontalis [55].

Detoxification genes number taken from studies of D. armandi[30], [100], [108].

Detoxification genes identified from genomes of D. ponderosae[52], A. glabripennis[53], and L. decemlineata[54].

Table 4

Detoxification genes up-regulated in D. valens and D. rhizophagus.

	Gen	Up-regulated	Log2FC	Log2CPM	p-value	FDR
D. valens
	CYP6DJ2	Control	1.06	6.17	7.01E-05	2.72E-03
	CYP334E1	With kairomones	−2.01	4.63	1.36E-06	1.15E-04
	CYP6DW23	With kairomones	−1.82	6.73	3.92E-13	3.32E-10
	CYP6DF1	With kairomones	−1.04	7.04	2.75E-03	4.60E-02
	CYP9Z20	With kairomones	−0.79	7.73	3.72E-04	9.89E-03
	CYP9Z52	With kairomones	−0.63	6.57	1.88E-03	3.51E-02
	CYP4BQ1	With kairomones	−1.02	4.76	2.63E-03	4.45E-02
	COEA3	With kairomones	−1.64	8.39	3.76E-12	2.27E-09
	COEA5	With kairomones	−0.80	5.64	5.70E-04	1.40E-02
	COEA11	With kairomones	−0.79	5.91	2.57E-03	4.40E-02
	COEA12	With kairomones	−1.17	10.22	1.99E-05	1.02E-03
	COEA15	With kairomones	−0.70	8.42	1.95E-03	3.59E-02
	COEA17	With kairomones	−1.67	5.29	8.07E-09	1.48E-06
	COEA18	With kairomones	−1.04	8.70	1.28E-06	1.10E-04
	COEA19	With kairomones	−0.92	6.30	3.88E-05	1.71E-03
	COEC1	With kairomones	−1.60	3.75	2.02E-03	3.68E-02
	COEC2	With kairomones	−0.62	8.39	1.74E-03	3.33E-02
	GSTo3	With kairomones	−1.38	6.50	2.22E-11	1.04E-08
	GSTs2	With kairomones	−1.23	6.12	5.91E-05	2.38E-03
	GSTs3	With kairomones	−1.18	7.38	1.75E-10	6.15E-08
	GSTe9	With kairomones	−2.06	3.23	1.81E-04	5.65E-03
D. rhizophaus
	CYP307A2	Control	1.07	7.77	3.41E-05	1.55E-03
	GSTz1	Control	1.38	5.79	2.18E-07	1.66E-05
	CYP6DE1	With kairomones	−2.10	7.83	1.62E-16	6.25E-14
	CYP6DE3	With kairomones	−1.92	9.49	5.83E-13	1.40E-10
	CYP6DJ1	With kairomones	−3.27	6.63	2.38E-30	2.01E-27
	CYP6DJ2	With kairomones	−2.31	4.43	1.20E-08	1.23E-06
	CYP9Z18	With kairomones	−1.00	8.34	9.27E-05	3.74E-03
	COEC2	With kairomones	−2.33	3.46	6.66E-06	3.59E-04
	GSTu2	With kairomones	−1.39	7.75	9.28E-08	7.44E-06

Cytochrome P450 monooxygenases

Transcriptome analysis of both species resulted in the identification of a total of 52 non-redundant P450 genes in D. valens (ACCN DvalP450s ON245918 - ON245969) and 51 in D. rhizophagus (ACCN DrhiP450s ON245867 - ON245917). The classification and naming of P450s was performed by the Nomenclature Committee (https://drnelson.uthsc.edu /CytochromeP450.html) (Table 3, Supplementary Table S5).

The complete DvalCYPs (41 genes) and DrhiCYPs (45 genes) varied from 446 to 576 amino acids, with molecular weights from 50.03 to 66.45 kDa, and pIs from 6.18 to 8.97 (Supplementary Table S6). The five conserved motifs present in insect P450s: the helix C motif (WxxxR), the helix I motif (Gx[ED]T[TS]), the helix K motif (ExLR), the PERF motif (PxxFxP[ED)RE), and the heme-binding motif (PFxxGxRxCx[GA]) [57] were identified in most of the P450s from both species (Supplementary Table S6).

In both species, the CYP2 and mitochondrial clans included the lowest number of P450s (Table 3). From the CYP2 clan, genes of the families CYP18, CYP303, CYP305, CYP306, and CYP307 were identified in both species, and only a gen of the subfamily CYP307B was recovered in D. rhizophagus. The mitochondrial clan was represented in both bark beetles by the families CYP49, CYP302, CYP314, CYP315, and CYP334; excepting CYP301, only recovered in D. valens. These families formed robust clusters in the phylogenetic tree and a close relationship with P450 genes from other beetles (Fig. 3).

Fig. 3

Maximum likelihood phylogenetic tree based of P450s from Dendroctonus valens (Dval) and Dendroctonus rhizophagus (Drhi). The analysis included P450s from bark beetles D. frontalis (Dfro), D. armandi (Darm) and D. ponderosae (Dpon), as well as those from Anoplophra glabripennis (Agla), and Leptinotarsa decemlineata (Ldec). The amino acid evolution model was LG + G (- lnL = 229583.22, G = 1.646, AIC = 463465.501). Branch support was calculated with the approximate likelihood ratio test (aLRT). Support values ≥ 80% are indicated on the branches with black point.

The CYP3 clan showed the greatest diversity with 27 genes in D. valens and 25 in D. rhizophagus, distributed in five families and 19 subfamilies (Supplementary Table S5). In both species, the CYP6 family was the most diverse within this clan, with a total of 16 genes in D. valens and 14 in D. rhizophagus. Only the families CYP6BS, CYP6CR, CYP6DE, and CYP6DH presented a different number of genes in D. valens (1,2,2,2) and D. rhizophagus (0,1,3,1). The phylogeny clustered P450s of these subfamilies with those from Dendroctonus species (Fig. 3), excepting DvalCYP6BS2 that was grouped with other beetle species. The second most diverse family was CYP9 with five genes in both species, distributed in three subfamilies (Supplementary Table S5). All CYP9 genes were integrated in exclusive clusters with P450s from other bark beetles in the phylogeny (Fig. 3). The families CYP345, CYP347, and CYP393 were recorded in both species and their members were clustered in the phylogeny with P450s from other beetle species (Fig. 3, Supplementary Table S5).

The CYP4 clan was represented by 14 genes in D. valens and 15 in D. rhizophagus, belonging to four families and 11 subfamilies (Table 3, Supplementary Table S5). The CYP4 family was the best represented within this clan, with 12 genes in D. valens and 13 in D. rhizophagus, distributed in the eight subfamilies; only the CYP4AA subfamily was exclusive from D. rhizophagus. The other families found were: CYP411, present in both species; CYP349, only recorded in D. valens, and CYP410, exclusive of D. rhizophagus.

Subfamilies CYP4AA, CYP4BR, CYP4G, CYP349B, CYP410A, and CYP411A from both bark beetle species clustered in the phylogeny with the P450s from other beetle species, but subfamilies CYP4BG, CYP4CV, CYP4BD, CYP4BQ, and CYP4EX were grouped exclusively with P450s from bark beetles (Fig. 3).

Carboxylesterases

A total of 37 non-redundant COE genes were identified in both species (ACCN DvalCOEs ON245830 – ON 245866, DrhiCOEs ON245793 – ON245829) (Table 3, Supplementary Table S5). The complete genes of both bark beetles (DvalCOEs 27 genes and DrhiCOEs 29 genes) varied from 428 to 1210 amino acids, with molecular weights from 58.96 to 137.55 kDa, and pIs from 4.14 to 8.09 (Supplementary Table S6). In most of the COEs from both species, the four conserved functional motifs present in insect COEs were identified (Supplementary Table S6). Among them are GxSxG, with the serine residue, x[DE]x, with the glutamic acid residue, and GxxHxx[DE], with the histidine residue; these three residues constitute the catalytic triad [58]. Additionally, the oxyanion hole domain HGGG was also identified [59].

The COEs from D. valens and D. rhizophagus clustered in the phylogeny within the three functional classes recognized by Oakeshott et al. [60]: 1) dietary and detoxification functions (27/28 genes), 2) hormone and semiochemical processing (5/4 genes), and 3) neurodevelopmental functions (5/5 genes), respectively (Fig. 4, Supplementary Table S5). Two big clades were integrated within the dietary and detoxification functions class: 1) The clade A (coleopteran xenobiotic metabolizing) represented by 21 DvalCOEs and 22 DrhiCOEs clustered mainly with bark beetle COEs; and 2) the clade C (microsomal and α-esterase) integrated by six COEs in both species, three of which were clustered with COEs from bark beetles, and other three with genes from other beetle species (Fig. 4). Despite the lack of variation in the catalytic triad residues, differences were observed in one or more residues within the four motifs. Such differences are reflected in the topology of the phylogeny and even in the specific subcellular location of some COEs in clade A (Supplementary Table S6).

Fig. 4

Maximum likelihood phylogenetic tree based of COEs from Dendroctonus valens (Dval) and Dendroctonus rhizophagus (Drhi). The analysis included COEs from bark beetles D. frontalis (Dfro), D. armandi (Darm) and D. ponderosae (Dpon), as well as those from Anoplophra glabripennis (Agla), and Leptinotarsa decemlineata (Ldec). The amino acid evolution model was VT + G (- lnL = 220061.043, G = 1.192, AIC = 443945.261). Branch support was calculated with the approximate likelihood ratio test (aLRT). Support values ≥ 80% are indicated on the branches with black point.

The hormone and semiochemical processing class comprises three clades: D (integument), E (β- and pheromone), and F (juvenile hormone). Within clade D, one COE was identified in both species, which form a single cluster with Dendroctonus species. Clade E was represented by three DvalCOEs and two DrhiCOEs, which were clustered with COEs from other beetles; whereas a single COE from clade F present in both species was clustered with COEs from other beetles studied (Fig. 4, Supplementary Table S5). The neurodevelopmental functions class was represented in both species by one gene in clades H, J, K, L, and M, which formed highly conserved clusters with other beetles in the phylogeny (Fig. 4, Supplementary Table S5).

Glutathione S-transferases

Twenty-eight glutathione S-transferase genes were identified in D. valens (ACCN: DvalGSTs ON245765 - ON245792), and 30 in D. rhizophagus (ACCN DrhiGSTs ON246919 - ON246941) (Table 3, Supplementary Table S5). The complete DvalGSTs (23 genes) and DrhiGSTs (26 genes) varied from 151 to 383 amino acids, with molecular weights from 22.61 to 43.68 kDa, and pIs varied from 4.44 to 8.68 (Supplementary Table S6). Three conserved motifs in all GSTs were identified: 1) the substrate selection motif (VPAL), 2) the GHS-binding motif (SNAIL/TRAIL), and 3) the substrate binding motif (GDxxxxAD) [61]. Additionally, we found other two motifs exclusive of the classes Delta and Epsilon (LYPx and RAxVxxRLxF) [62] (Supplementary Table S6).

The phylogenetic analysis showed that 25 DvalGSTs and 27 DrhiGSTs were distributed in the seven recognized cytosolic GSTs classes (Theta, Zeta, Omega, Sigma, Delta, and Epsilon). From Zeta and Theta classes was found one GST in each species, which were clustered with members of these classes from other beetles in the phylogeny. Three Omega GSTs were found in both species, two of which were clustered with the genes of bark beetles; the other one was clustered with GSTs from other beetles. From Sigma class, six DvalGSTs and seven DrhiGSTs were found, which integrated exclusive clusters with genes from other Dendroctonus species. Two Delta GSTs in D. valens and three in D. rhizophagus were identified, two of them, common for both species, were clustered with GSTs from other beetles, the other DrhiGSTd3 was associated with GSTs from other beetles. Twelve Epsilon GSTs were recorded in both species, most of them were clustered with those from Dendroctonus species (Fig. 5).

Fig. 5

Maximum likelihood phylogenetic tree based of GSTs from Dendroctonus valens (Dval) and Dendroctonus rhizophagus (Drhi). The analysis included GSTs from bark beetles D. frontalis (Dfro), D. armandi (Darm) and D. ponderosae (Dpon), as well as those from Anoplophra glabripennis (Agla), and Leptinotarsa decemlineata (Ldec). The amino acid evolution model was VT + G + I (- lnL = 220061.043, G = 1.800, I = 0.003, AIC = 41038.861). Branch support was calculated with the approximate likelihood ratio test (aLRT). Support values ≥ 80% are indicated on the branches with black point.

The remaining sequences of both species were associated with unclassified GSTs documented in other beetles (Table 3; Fig. 5). Two of them were related to GSTs from Delta class reported in other beetles. The other unclassified GSTs were clustered with the GSTs of Dendroctonus species, which are closely related to the Omega GSTs (Fig. 5).

Intra-and interspecific comparison of detoxification gene expression

A total of 21 detoxification genes were differentially expressed in D. valens and only nine in D. rhizophagus. Twenty (6 P450s, 10 COEs and 4 GSTs) and seven (5 P450s, 1 COE and 1 GST) of them were up-regulated in response to kairomones, respectively. The up-regulated genes in D. rhizophagus had a higher Log2FC than the up-regulated genes in D. valens, even though the number of up-regulated genes in D. valens is greater than in D. rhizophagus (Table 4).

A comparison of the expression profile between the detoxification genes recorded in D. valens (1 1 7) and D. rhizophagus (1 1 8) revealed that only 77.78% (34 P450s, 32 COE and 25 GST) and 76.27 % (35 P450, 31 COE and 24 GST) had detectable TPM values (>5 TPM), respectively, in both conditions. In response to the kairomone blend, DvalCOEA20 (5.44) and DrhiCOEF1 (6.09) had the lowest TPM value, while DvalCOEA12 (1099.56) and DrhiCOEA12 (1048.40) had the highest value. Interestingly, most of the genes presented a species-specific transcriptional response (Fig. 6, Supplementary Table S6).

Fig. 6

Expression profile of detoxification genes from D. valens and D. rhizophagus based on TPM values. Genes expression levels were indicated with scale of color using a mean TPM value of two replicates.

Particularly, in both species, genes of the CYP4 and CYP6 families, from P450s, as well as of the clade A, from COEs, and of the Sigma class, from GSTs, displayed higher TPM than other genes. However, most of these did not present significant differences with respect to the control samples (Fig. 6, Supplementary Table S6).

Discussion

Comparative genes expression profiles between D. valens and D. rhizophagus

Our findings support a reliable gut transcriptome assembly in both bark beetles, whose transcripts encode full-length protein sequences. Both bark beetles showed difference in the number of up-regulated unigenes, being major in D. valens (739) than D. rhizophagus (322), in response to the kairomones blend (α-pinene, β-pinene and 3-carene). These results suggest that these bark beetles have a different adaptive physiological response to the harmful effects of the kairomones/terpenes, as a consequence related to their life histories and ecological strategies. D. valens is a gregarious, multivoltine, and non-aggressive species in its original distribution range in North America, where it colonizes weakened or dying adult trees or stumps from around 40 pine species [63]; in addition, it does not produce massive attacks despite several couples arriving at the same tree. In contrast, D. rhizophagus is a non-gregarious and univoltine species, that does not produce massive attacks, yet it is an aggressive species endemic to the Sierra Madre Occidental in Mexico where only one couple, occasionally two, colonize and kill young hosts (seedlings ≥ 1.5 cm diam. base to saplings < 8 cm diam. at 1.4 m height and <3 m tall) of 11 pine species [64], [65].

The major transcriptional response recorded in D. valens compared to D. rhizophagus might also be indicative of a major tolerance to a wide range of chemical compounds produced by adult trees, whose terpenes composition and concentration varies considerably among tree species, age, season, and geographic space [6]. In contrast, the transcriptional response observed in D. rhizophagus might be associated to a major efficiency and specificity toward chemical compounds present in the seedlings of trees.

Our findings agree with other studies that compare the transcriptional response of phytophagous insects classified as generalist and specialist species, based on the diversity and number of host plants they use [66]. These studies have documented that generalist species (those that commonly use numerous hosts of the same or different genus or even different families) have higher transcriptional response than specialist species (those that use a few plants, regularly of the same genus) to overcome the defensive chemical compounds of host plants. In contrast, specialist species possess more specific and efficient enzymatic systems. In addition, substantial changes in the expression of genes, directly or indirectly related to xenobiotic detoxification, have been observed between both behaviors [67], [68], [69].

It is expected that genes of the three phases of the detoxication process (P450s, COEs, GSTs and ABC transporters) involved directly or indirectly in the metabolism and excretion of terpenes in D. valens and D. rhizophagus are those that show an increment in their transcriptional response. Yet, it is notable that up-regulated unigenes in D. rhizophagus (P450s and some COEs) were different to those in D. valens (COEs, GSTs and ABC transporters). These results indicate that these species have specific capacities to metabolize kairomones.

Interestingly, phase I and II reactions as well as their combination, can participate in the bioactivation of xenobiotics leading to the production of chemically reactive intermediates that bind covalently to DNA or tissue proteins producing adducts potentially toxic to the cells and organs of the insect [70]. In addition, other molecules resulting of the detoxification metabolism might produce direct cellular damage, jeopardizing the integrity of cells and organs [71]. This suggests that terpenes metabolism comes accompanied of integral hormetic responses (adaptive responses that increase the resistance of the cells and/or organism to severe and/or lethal stress) that result of compensatory mechanisms, such as different metabolic nertworks (KEEG pathways) and/or functional categories (GO terms) recorded in these bark beetles (Fig. 2).

Studies carried out in phytophagous insects suggest that hormetic responses triggered by xenobiotics metabolism induce metabolic networks that comprising enzymes of the phase I (COEs), phase II (GST), and antioxidant, as well as cellular repair mechanisms and nutrients acquisition that ensure the supply of energy and the survival of the organism [72], [68], [69].

The diversity of P450 genes found in D. valens (52) and D. rhizophagus (51) is similar to that found in the transcriptome of D. armandi (46) but lower than those in the genomes of D. ponderosae (76), L. decemlineata (90), and A. glabripennis (128) (Table 3).

The genes number in the families of the mitochondrial and CYP2 clans was minimal in both bark beetles (Table 3), showing a low nucleotide divergence between them (Fig. 3). This is observed in the phylogeny, where these genes are clustered according to essential physiological processes where they participate, such as development and reproduction [57], molting and metamorphosis [73], and xenobiotics metabolism [74].

Among the genes of mitochondrial and CYP2 clans highlight the Halloween genes: CYP302A1 (Disembodied or Dib), CYP306A1 (Phantom or Phm), CYP307A2 (Spookier or Spok), CYP314A1 (Shade or Shd) and CYP315A1 (Shadow or Sad), as well as CYP305F1 and CYP18A1, all of which are involved in the synthesis and regulation of ecdysteroids [75], [73], [76]. The remaining P450s found in this study have been associated to diverse functions in other insects, such as the regulation of cuticular hydrocarbons (CYP49A1), cuticle formation (CYP301A1), signal metabolism in sensorial organs (CYP303A1) and gonad maturation (CYP307B1, Spookiest or Spot) [77], [78], [79], [80]. The differences in expression levels observed in orthologous genes of these bark beetles, especially DvalCYP49A1, DvalCYP334E1 and DrhiCYP307A2 (Fig. 5) could be result of the physiological state in which the insects of both species were collected.

The number of genes of CYP3 clan recorded in D. valens (27) and D. rhizophagus (25) was practically the same; however, they had important differences in expression profiles. While the functional role of the P450 enzymes from D. valens and D. rhizophagus is not yet known, specific studies have analyzed and recorded significant differential expression in specific P450s of these families in several Dendroctonus species [20], [22], [30]. For example, in D. rhizophagus and D. valens, high expression levels have been demonstrated in CYP6DJ1, CYP6DJ2, CYP6BW5 and CYP6DG1 after insects were stimulated with host terpenes vapors [40], [41], as well as during different colonization phases in the field [19] and during the early hours of drilling and settling in the host tree by D. rhizophagus [21]. The up-regulation of DrhiCYP6DJ1, DrhiCYP6DJ2, DrhiCYP6DE1 and DrhiCYP6DE3 observed in this study, suggests that these enzymes might be directly associated to the detoxification process.

In D. ponderosae, it has been demonstrated that CYP6DJ1 metabolizes terpinolene and limonene to alcohols and an epoxide [25]; CYP6DE1 metabolizes (±)-α-pinene, (±)-β-pinene and (+)-3-carene to different oxygenated monoterpenes, being the trans-verbenol the main product derived of the activity of this enzyme with both enantiomers of α-pinene [24], whereas CYP6DE3 produces trans-verbenol only from (+)-α-pinene [32]. Interestingly, in D. rhizophagus, trans-verbenol has been identified as a probable aggregation pheromone [37], which might explain the up-regulation of genes DrhiCYP6DE1 and DrhiCYP6DE3 in this species.

On the other hand, the up-regulation of DvalCYP6BW23 may be related to the oxidization of various diterpene acids as observed in CYP6BW1 and CYP6BW3 from D. ponderosae [25]. On the other hand, it has been demonstrated that D. valens produces frontalin (sexual and aggregation pheromone) and exo-brevicomin (anti-aggregation pheromone) [38], [39], DvalCYP6CR1 and DvalCYP6DH2 were not up-regulated, as their orthologs in D. ponderosae, which are involved in the synthesis of both pheromones. The first one, epoxidizes (Z)-6-nonen-2-one to 6,7-epoxynonan-2-one, a precursor of the pheromone exo-brevicomin [18], [81], [32], [82], and the second is considered an intermediary candidate in the frontalin biosynthesis [32]. Another gene probably involved in the exo-brevicomin biosynthesis in D. ponderosae [83], [32] is CYP6DF1, whose ortholog DvalCYP6DF1 was up-regulatated in this study.

It should be noted that none of these orthologs are present in D. rhizophagus; their absence might explain why studies on chemical ecology and expression of mevalonate pathway genes have not recorded frontalin or exo-brevicomin in this species [84]. Duplicates of these genes (CYP6CR2 and CYP6DH3) are present in D. valens, D. rhizophagus, and D. ponderosae, but there is no evidence that they are involved in the synthesis of any semiochemical.

Another family involved in the metabolism of xenobiotics is the CYP9 [57]. In this study, genes of the subfamilies CYP9AP, CYP9AZ, and some members from CYP9Z were identified in both bark beetle species. The up-regulation DvalCYP9Z20, DvalCYP9Z52 and DrhiCYP9Z18 could be the result of the stimulus of the different kairomones, because experimental studies have recorded the significant induction of genes CYP9Z18 and CYP9Z20 in both bark beetles and D. armandi with different terpenes and under different experimental conditions [40], [41], [30]. Duplicates of these genes (CYP9T2 and CYP9T3) have also been reported in the Ips-bark beetles, which participate in the de novo synthesis of the aggregation pheromones ipsdienol and ipsenol, via the mevalonate pathway from myrcene, and whose capacity to transform other terpenes has also been demonstrated [85].

The presence of other less diverse families of the CYP3 clan (CYP345, CYP347, and CYP393) also suggests its involvement in detoxification processes. It has been shown that the CYP345E2 from D. ponderosae oxidizes (±)-α-pinene, (±)-β-pinene, (+)-3-carene, terpinolene, (±)-limonene, and (-)-camphene [23], and CYP345F1 has been reported as a probable candidate involved in some steps of the biosynthetic pathways of frontalin, through co-expression networks [32]. These families are also present in other insects, such as T. castaneum, whose expression have been linked to the xenobiotic’s metabolism [86].

We hypothesized that the main function of CYP6, CYP9 CYP345, CYP347, and CYP393 families in D. valens and D. rhizophagus is related to terpenes detoxification. Thus, gene duplication in these families might be primarily linked with a subfunctionalization process to increase the detoxification capacities of these bark beetles, and, secondarily, participate in compensatory mechanisms or pheromone biosynthesis through different metabolic pathways, as it is evident in other bark beetle species [87].

With relation to CYP4 clan a similar number of transcripts was also recorded in D. valens (14) and D. rhizophagus (15). This clan comprises subfamilies of the highly divergent CYP4 family [57]. In fact, members of the CYP4 family are involved in several physiological functions in phytophagous insects, such as endogenous compounds biosynthesis [88], odorant degradation [89], fatty acids metabolism and cuticular hydrocarbons synthesis [90], as well as cuticle formation and pigmentation [91]. Members of this family are integrated in well-defined groups in the phylogeny, and their involvement in conserved physiological functions is suggested (Fig. 3).

Despite that the physiological function of many CYP4 in bark beetles is not yet known, several studies have expressed genes of different subfamilies (e.g., CYP4C, CYP4G, CYP4EX) in different tissues and systems of the genus Dendroctonus under different experimental conditions [40], [41], [30]. In fact, it has also been suggested that the genes CYP4CV2 and CYP4EX1 from D. ponderosae are also probably involved in exo-brevicomin biosynthesis [32]. Interestingly, our results revealed the presence of two genes from the subfamily CYP4C, CYP4CV1 and CYP4CV2, both in D. valens and D. rhizophagus, which are present in other Dendroctonus species, as well as in the coffee berry borer bark beetle, Hypothenemus hampei [27], [30], [92]. Studies in cockroach strongly suggest that CYP4C7 is involved in the biosynthesis of juvenile hormone (JH) through the metabolism of JH-precursors and JH itself within the corpora allata (CA) during specific developmental stages [88]. In bark beetles, it has been demonstrated that phloem feeding during trees colonization triggers different physiological events, such as JH-biosynthesis in the CA, flight muscle degeneration, and the de novo pheromone biosynthesis via the mevalonate pathway which is regulated by JH [93]. By analogy, we hypothesized that enzymes CYP4CV1 and CYP4CV2 together with JH-esterases could regulate the levels of this hormone in these Dendroctonus species.

The record of CYP4G55, CYP4G56 and CYP4G110, in both D. valens and D. rhizophagus, is another noteworthy result of this study. These three genes have been reported in D. armandi [30], but only CYP4G55 and CYP4G56 genes in D. ponderosae [27]; however, a targeted search in the genome of this last species allowed us to identify a fragment of the CYP4G110 gene (ERL89343.1, 509aa). Orthologs of these three genes have also been found in H. hampei (CYP4G153, CYP4G155, CYP4G154) [90]. Based on the phylogenetic relationship and in silico characterization of the enzymes encoded by these genes (data not shown), we hypothesize that a double gene duplication occurred in Dendroctonus species. The first one between CYP4G55/56 and the second between CYP4G55/CYP4G110; a similar event has been suggested for CYP4G155/154/153 in H. hampei [90].

This double duplication is relevant for Dendroctonus species and other bark beetles, because these genes participate in the terminal steps of cuticular hydrocarbon biosynthesis, by oxidizing long chain primary alcohols to aldehydes, and subsequently decarbonylating aldehydes to hydrocarbons and carbon dioxide [83], [94]. In addition, cuticular hydrocarbons are highly specific and might function as cues for gender recognition [95]. They also, together with the thickness and the degree of cuticular sclerotization, are a determining physical barrier to avoid dehydration [96] and penetration of toxic terpenes from the host tree into the insect’s body or within the gut cells, where the presence of sclerotized plaques and grooves have been documented [16], allowing the enzymes involved in the detoxification process to adjust and increase their metabolic capacity [97]. Additional studies on CYP4G110 gene could confirm its participation in this metabolic process.

It is noteworthy that CYP4AA1 had been recovered only in D. rhizophagus. The CYP4A subfamily members include enzymes hydroxylating fatty acids at the terminal ω-position, an unusual hydroxylation of medium-chain fatty acid, in contrast to the P450 enzymes that have preference for non-terminal hydroxylation and⁄or epoxidation. The CYP4AA1 acts as a fatty acid ω-hydroxylase in pheromone biosynthesis in honeybee which regulates caste and social condition [98]. However, in several beetle species where orthologs of this gene have been reported, its function remains unclear.

Lastly, although the functions of CYP349B1, CYP410E1, and CYP411A1 genes are unknown in bark beetles, their presence in Dendroctonus spp. suggests, that they also play an important role in the detoxification process of these species, because it has been shown in D. armandi that these genes have differential expression in different development stages under stimuli with terpenes [30].

Despite their potential relevance in detoxification processes, the functional characterization of COEs in bark beetles and specifically in Dendroctonus species has been poorly studied. Our findings show a diversity of COEs in D. valens (37) and D. rhizophagus (37) comparable with that identified in the genome of D. ponderosae (51), but lower than reported in the genome of L. decemlineata (80) and A. glabripennis (107) (Table 3).

The dietary and detoxification functions class is widely represented in all insects and include members of the clades A and C. Particularly, the COEs of clade A constitute a diverse group involved mainly in the xenobiotic metabolism, and some of them in the metabolism of insecticides in several species such as L. decemlineata [99]. Members of this clade integrate, exclusive and well-defined clusters with COEs of bark beetles in the phylogeny. The formation of these cluster is due to specific amino acids changes in motifs of these enzymes, which are synonymous changes that not affect the catalytic triad (SEH). The changes in Dendroctonus species occurs mainly in the GxSxG motif, allowing the 3D structure of these proteins to produce the nucleophilic elbow (an active catalytic site able to interact with different substrates) [59].

The up-regulation of 10 COEs in D. valens after exposure to the kairomones-blend, suggests that these COEs have played a fundamental role in the metabolism of terpenes; however, these genes in D. rhizophagus are not up-regulated. In D. armandi and D. ponderosae, it has been shown the up-regulation of some of these COEs after the insects were fed on phloem or stimulated with monoterpene vapors [29], [100], [22].

In the case of COEs of clade C, despite their function has remained unclear, the formation of exclusive clusters of orthologs in the Dendroctonus species analyzed, showing a low nucleotide divergence between them, suggesting an important role in detoxification and digestive processes. In fact, DvalCOEC1, DvalCOEC6, DrhiCOEC2 were up-regulated in presence of the blend of kairomones.

Regarding the hormone and semiochemical processing class, integrated by clades D, E, and F, it has been suggested that members of these clades participate in the inactivation of pheromones or kairomones during signaling, influencing the reproductive behavior of insects [101], [102], [103]; nonetheless, there is no evidence that these enzymes participate in semiochemicals processing in Dendroctonus species.

The members of the clade D are known as integument esterases implicated in the processing of olfactory signals (pheromone and kairomone) [104], which may explain why, no member of this clade was up-regulated in both species. Regarding the members of clade E, these have been associated with a variety of functions ranging from the processing of sex pheromones and organophosphates (OPs) resistance [63], [103], [104]. In this study, DvalCOEE2 and DrhiCOEE2 did not show differential expression, but its ortholog in D. armandi (DaCarE1) showed up-regulation in females exposed to (+)-3-carene [100].

On the other hand, it is known the role of COEs of clade F in the processing of JH [99]. In fact, DvalCOEF1 and DrhiCOEF1 correspond to the juvenile hormone esterase deposited in the GenBank. The same gene of this esterase, DaCarE4, has been identified in D. armandi, which was significantly up-regulated in both sexes after feeding and exposure to (+)-3-carene and turpentine at 8 and 24 h [100], [22].

Finally, the neurodevelopmental functions class that include members of the clades H, J, K, L and M showed to be phylogenetically conserved, because they were integrated with proteins of the species of beetles analyzed. These genes are mostly catalytically inactive, because some studies performed in D. melanogaster and A. mellifera have shown that they act as specific cell adhesion molecules during nervous system development [60], [105].

The number of cytosolic GSTs found in D. valens (28) and D. rhizophagus (30) was similar to those identified in the genome of other beetles, such as D. ponderosae (32), L. decemlineata (32), and A. glabripennis (37) (Table 3). Yet, important differences in the Sigma, Delta, and Epsilon classes were observed (Table 3). The GSTs Zeta and Theta classes from D. valens and D. rhizophagus are orthologous to those of other beetle species analyzed. The Zeta GSTs have been related to the degradation of tyrosine and phenylalanine once the cuticle sclerotization process is completed [106]. By analogy, the up-regulation of DrhiGSTz1 in sample control might be associated with the maturation of the adult cuticle in D. rhizophagus. Some reports also associate the activity of GSTz to the detoxification of chloride-containing xenobiotics [107].

The presence of three Omega GSTs in Dendroctonus species along with the specific up-regulation of DvalGSTo3, suggest that these GSTs might play an important role in the adaptation of these species to their toxic environment. The ortholog of this gene in D. armandi (DarmGSTo2) was significantly up-regulated in both sexes with (−)-α-pinene and (−)-β-pinene [108]. Members of this class are involved in protection against oxidative stress in Apis cerana and Rhopalosiphum padi. In this last species, it has observed that the expression of GSTO1 occurs when insects are exposed to different insecticides [109], [110].

In a similar way, the presence of Sigma GSTs in D. valens (6) and D. rhizophagus (7) also suggests a possible role in the adaptation of these bark beetles to toxic subcortical environments. It has been proposed that GSTs of this class participate in protection against oxidative stress by exhibiting conjugation activity for lipid peroxidation products and carbaryl insecticides [111], [112], [113]. This function has been related to the expansion of this class in Nasonia vitripennis, where it has been hypothesized that these enzymes can protect this species against the reactive oxygen molecules produced by the progressive death of its hosts [104].

Given that bark beetles lead their host trees until death, the Sigma GSTs of bark beetles may play a similar functional role in protection against oxidative stress involved in the detoxification process. In fact, DvalGSTs2 and DvalGSTs3 were up-regulated in this study, and analyses of RNA-seq and proteome profiles from other Dendroctonus species, fed on fresh phloem and exposed to different terpenoids, showed a significant increment in the expression of members of Sigma GSTs class [29], [114]. In addition, Sigma GSTs in D. armandi have showed variation in their expression under different conditions [108], [115], [22], which supports that these enzymes play an important role in reducing the negative impact of terpenoids on these beetles.

With respect to the Delta, Epsilon and Unclassified classes, their members have been mainly associated with resistance to insecticides and secondary metabolites of plants [116]. In our comparative analysis with these two bark beetles, members of these three classes are represented in more than 50% of the whole GST superfamily. Genome studies in coleopterans have documented multiple isoforms derived from alternative splicing [14], which is explained as an adaptive response of coleopterans to selective pressures of their environment. In these classes, two extra motifs have been identified, which suggested they are binding sites to substrates that allow their interaction with a great variety of compounds [65].

Our phylogeny shows that the Delta class might have the same function in beetles, because clusters were integrated independently of beetle species. On the other hand, the enzymes of the Epsilon class apparently have an exclusive role in Dendroctonus spp., because they form clusters independent from other beetles. The integration of these clusters in Dendroctonus bark beetles could be result of specific selective pressures in its subcortical habitat. Important variations in the expression of Delta and Epsilon GSTs in D. armandi have been observed between sexes, development stages, and treatments (e.g. feeding, terpenoid and combined treatments) [108], [115], [22]. It is noteworthy that none of these genes was up-regulated in D. valens and D. rhizophagus.

However, DrhiGSTu2 gene, belonging to the unclassified class, was up-regulated in the presence of kairomones blend. In Locusta migratoria a considerable increment in mortality occurred in nymphs treated with carbaryl and chlorpyrifos by silencing LmGSTu1, suggesting that this gene plays a significant role in the detoxification of these insecticides [117]. In Bombix mori, bmGSTu2 conjugates glutathione to 1-chloro-2,4-binitrobenzene, as well as metabolizes diazinon, one of the most used organophosphate insecticides, and it is closely related to epsilon-class GST [118]. In L. decemlineata, LdGSTu1 could catalyze the conjugation of GSH to both CDNB (1-chloro-2, 4-dinitrobenzene) and PNA (p-nitrophenyl acetate), and inhibition assays demonstrated that the enzymatic conjugation of GSH to CDNB was inhibited by multiple pesticides suggesting a potential function of LdGSTu1 in xenobiotic adaptation [119].

Conclusion

Our results showed that the expression profile of P450s, COEs, and GSTs genes in two sibling species of bark beetles exposed to the same kairomones blend is species-specific. In addition, the stimuli with this blend induce the expression of genes related to different metabolic pathways and functional terms, which are also different between these bark beetles. These metabolic and functional processes show that the adaptive responses of these bark beetles to terpenes are complex and comes accompanied by an integral hormetic response. Further studies are needed to depth in the knowledge of the differentially expressed genes in these beetles and to get a better understanding of their integral adaptive response.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.