Comparative transcriptomics of mountain pine beetle pheromone-biosynthetic tissues and functional analysis of CYP6DE3
© The Author(s). 2017
Received: 30 November 2016
Accepted: 11 April 2017
Published: 20 April 2017
The mountain pine beetle (MPB, Dendroctonus ponderosae Hopkins) is a highly destructive pest of pine forests in western North America. During flight to a new host tree and initiation of feeding, mountain pine beetles release aggregation pheromones. The biosynthetic pathways of these pheromones are sex-specific and localized in the midgut and fat body, but the enzymes involved have not all been identified or characterized.
We used a comparative RNA-Seq analysis between fed and unfed male and female MPB midguts and fat bodies to identify candidate genes involved in pheromone biosynthesis. The 13,407 potentially unique transcripts showed clear separation based on feeding state and gender. Gene co-expression network construction and examination using petal identified gene groups that were tightly connected. This, as well as other co-expression and gene ontology analyses, identified all four known pheromone biosynthetic genes, confirmed the tentative identification of four others from a previous study, and suggested nine novel candidates. One cytochrome P450 monooxygenase, CYP6DE3, identified as a possible exo-brevicomin-biosynthetic enzyme in this study, was functionally characterized and likely is involved in resin detoxification rather than pheromone biosynthesis.
Our analysis supported previously characterized pheromone-biosynthetic genes involved in exo-brevicomin and frontalin biosynthesis and identified a number of candidate cytochrome P450 monooxygenases and a putative cyclase for further studies. Functional analyses of CYP6DE3 suggest its role in resin detoxification and underscore the limitation of using high-throughput data to tentatively identify candidate genes. Further functional analyses of candidate genes found in this study should lead to the full characterization of MPB pheromone biosynthetic pathways and the identification of molecular targets for possible pest management strategies.
KeywordsTranscriptomics RNA-Seq Bark beetle Dendroctonus Pheromone P450
Various high-throughput studies, including a draft MPB genome , and various transcriptomic [7, 8] and proteomic  analyses have been directed at this highly destructive coniferous forest pest in order to better understand its physiology and to assist development of new management strategies. Aggregation and anti-aggregation pheromones in the MPB are synthesized in the midguts and fat bodies in a sex- and feeding status-specific manner . Therefore, mRNAs encoding enzymes involved in pheromone-biosynthetic pathways may also be differentially expressed based on sex, feeding status, and tissue. Recently, Keeling et al.  reported transcriptome, proteome, and metabolome-level responses of unfed MPB treated with juvenile hormone (JH) III. JH III stimulates aggregation pheromone component production in bark beetles [5, 12–15], including frontalin and trans-verbenol in MPB, and several “pheromone-biosynthetic gene” candidates were identified by this study. However, the dynamics of the transcriptional response of JH III-treated insects can differ from those of fed insects [15–17] and the study noted that JH III treatment did not affect exo-brevicomin production. Thus, JH treatment may not induce differential expression of all genes involved in pheromone biosynthesis for the MPB.
Here, we report the results of an RNA-Seq study of MPB midguts and fat bodies based on sex and feeding status. We employed network, differential gene, and gene ontology enrichment analyses to isolate candidates involved in pheromone biosynthesis. Additionally, functional analysis of a candidate gene, CYP6DE3, identified by our RNA-Seq analysis for exo-brevicomin biosynthesis, indicates a role in resin detoxification rather than pheromone biosynthesis, reinforcing the need for functional assays to confirm possible roles of enzymes identified via bioinformatics approaches.
Sections of mountain pine beetle-attacked lodgepole pine, Pinus contorta, were collected from Truckee, CA, USA, near the Northwoods Clubhouse in the Tahoe Donner housing subdivision (approx. N 39°20′37″ W 120°12′54″) on September 30, 2013. The beetles overwintered in the bolts and emerging adult beetles were collected and sexed as reported previously [3, 18]. For feeding experiments, fresh lodgepole pine bolts were obtained from the Whittell Forest, located in the Carson Range on the east slope of the Sierra Nevada (approx. N 39°16′29″ W 119°52′43″) in June, 2014 and stored at 4 °C prior to use. Females were fed by drilling small holes just beneath the bark, inserting the beetles head first, stapling a wire mesh over the occupied hole, incubating the bolt vertically for 24 h in the dark at room temperature and collecting the live beetles by stripping the bark. Fresh frass indicated that the beetles had fed. Males were fed using the same method except females were first placed head first into the holes for 24 h, followed by insertion of the males head first for 24 h and subsequent collection and sexing of the live beetles. Unfed beetles were incubated in 2 oz plastic cups with perforated lids in a dark drawer kept humid with small flasks filled with water and a paper towel for 24 h. All beetles were dissected following treatment to collect midgut and fat body tissue. Four replicates of pooled tissue from 10 beetles were collected for each of the four treatments (fed and unfed males and females) for a total of 16 samples. Midguts and fat bodies were immediately frozen in liquid nitrogen and stored at -80 °C for subsequent RNA extraction.
RNA extraction and sequencing
Total RNA was isolated from the midguts and fat bodies pooled from 10 beetles per sample using an RNeasy Plant Mini Kit from Qiagen (Valencia, CA) and treated with RNase-Free DNase Set from Qiagen as described by the supplier’s manual. RNA was quantified using Quant-iT™ RiboGreen® reagent and a Labsystems Fluoroskan Ascent fluorescence plate reader. RNA integrity for each sample was determined using an Agilent 2100 Bioanalyzer and a Eukaryote Total RNA Nano Series II chip. Only samples with an RNA Integrity Number of eight or higher were used for sequencing . RNA was precipitated and provided to the Georgia Genomics Facility (GGF) for library preparation and sequencing. GGF confirmed the quality of the total RNA using an Agilent 2100 Bioanalyzer, prepared barcoded cDNA libraries using a Kapa Stranded mRNA-Seq Kit (Wilmington, MA), and sequenced them on the Illumina NextSeq 500 using paired-end sequencing with a NextSeq 2x75 High Output Flow Cell.
Sequence quality control
Sequence quality for each sample was characterized using FastQC (v. 0.11.2; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Sequence pairs were trimmed and filtered for nucleotide-base quality and Illumina adapter sequences using Cutadapt v. 1.8.1 , with options set as follows: trimming of low-quality (phred quality ≤ 10) and “N” base calls from both ends of each read and removal of sequences with trimmed length < 35 nt.
Sequence alignment and expression quantification
Before assembly, trimmed sequence pairs were compared to one another using the MaSuRCa ‘superreads.pl’ script, and if found to intersect (minimum k-mer 41), were combined into single “super-reads” . Read pairs with no intersection were retained as separate paired-end reads. Single reads and read pairs were aligned to the Ensembl Metazoa (release-25; ) D. ponderosae reference genome  using the HISAT spliced read alignment tool (v. 0.1.6-beta; ). The coordinates of each known Dendroctonus gene and its exons were extracted from the Ensembl Metazoa Gene Transfer Format (GTF) file and supplied to HISAT at time of alignment (via HISAT’s ‘--known-splicesite-infile’ option). Resultant alignments were compressed from the sequence alignment/map (SAM) format to the binary BAM format . Upon alignment, the raw counts of reads and read pairs aligned to each gene were totaled using the featureCounts tool of the subread package (v. 1.4.6-p1; ). Reads were counted once per pair and summarized for gene loci, with only read pairs aligned to a unique transcribed location included in the count totals.
Transcripts underwent two filtering steps. First, those with no counts for thirteen or more out of the sixteen (75%) samples were excluded. Then, all transcripts with less than 10 fragments (counts) observed in all sixteen samples were excluded. Data were normalized using the standard median ratio method for RNA-Seq data . Principal component analysis (PCA) was performed on the normalized and filtered zero-centered counts per million data using singular value decomposition to validate clear separation between gender and feeding status of the biological replicates of MPBs.
Dendroctonus ponderosae gene descriptors and annotated Interpro protein domains  were obtained from Ensembl Metazoa, via the BioMart interface . Interpro2GO file (v. 2016/03/19 11:04:26) was used to map Interpro IDs to gene ontology (GO) terms.
Validation of RNA-Seq data by quantitative PCR
To validate RNA-Seq data, the transcript levels of 15 genes (Table 2) were examined by quantitative reverse transcriptase PCR (qRT-PCR). Genes were chosen by their notable differential expression between feeding states in male beetles. Aliquots consisting of approximately 500 ng of total RNA from a subset of samples (fed male replicate 4, unfed male replicate 4, fed female replicate 4, unfed female replicate 4, fed male replicate 1, unfed male replicate 3, unfed female replicate 2) were used to make cDNA using iScript™ Reverse Transcription Supermix (Bio Rad, Hercules, CA). PCR was conducted in a 20 μL reaction consisting of iTaq™ Universal SYBR® Green Supermix (Bio Rad) and 2 μL of template for 40 cycles of 95 °C for 5 s and 60 °C for 30 s on a Bio Rad CFX96 Real-Time PCR Machine (Bio Rad). Primers designed to amplify specific transcripts of these genes were designed using IDT Primer Quest and melt curves were produced to ensure primer specificity and proper PCR temperature cycling parameters (Additional file 1). For each cDNA sample the PCR reactions were conducted in triplicate and relative target gene expression was normalized to that of YQE_05788, which encodes ribosomal protein S3P. Ribosomal protein S3 is a relatively more stable normalizing gene for qRT-PCR in another beetle, Tribolium castaneum, compared to the more routinely used actin or tubulin genes . Fold change was calculated for each normalized gene in relation to the expression of the unfed male treatment using the 2-ΔΔCT method . For each gene, Pearson and Spearman Correlation Coefficients were computed between the seven samples measured by qRT-PCR and RNA-Seq.
Gene co-expression networks are node-edge graphs. Nodes represent genes that are connected by edges if there is an association between genes as defined by a co-expression measure . Structural components of co-expression networks are used to identify densely connected subgraph, called gene modules. Genes within a module share similar expression patterns, thus they are hypothesized to have similar gene function, to share pathway membership, or to be co-regulated. A co-expression network of the filtered and normalized counts per million of 11,342 mountain pine beetle genes was generated via petal, a co-expression network construction and analysis tool . The entire dataset of 11,342 genes over 16 measures was loaded, along with a list of previously confirmed and hypothesized pheromone biosynthetic genes. No other input was specified.
Differential gene analysis
Differential gene expression between the feeding conditions and the genders were examined using DESeq2 . Four comparisons, male fed (MF) vs. female fed (FF), MF vs. male unfed (MU), MU vs. female unfed (FU) and FF vs. FU, were considered using simple contrasts. A multiple testing correction was performed for each of the four comparisons to adjust for the false discovery rate . The two other possible comparisons (MU vs. FF and MF vs. FU) were not considered because they are less likely to inform regarding putative pheromone-biosynthetic genes. Genes with absolute value of the log2 fold change greater than one and an adjusted p-value less than 0.01 were retained for further analysis. Venn diagrams were prepared within R to visualize the intersection of the statistically-significant differentially-expressed genes between the considered comparisons (Additional file 2).
Gene ontology enrichment analysis
Gene Ontology (GO) enrichment analyses of statistically-significant differentially-expressed gene groups were conducted to identify over-represented molecular functions and metabolic processes. BiNGO (v. 3.0.3)  within Cytoscape (v. 3.3.0) and GO file (v. 1.2 2016/03/01) were utilized. GO terms with adjusted significance values less than 0.05 upon a Benjamini-Hochberg adjustment  were considered for further investigation.
Recombinant CYP6DE3 expression
RNA was extracted from two whole beetles using a RNeasy Plant Mini Kit (Qiagen) as per the manufacturer’s instructions. First strand cDNA was synthesized using Superscript III Reverse Transcriptase (Invitrogen, Carlsbad, CA) as per the manual. The CYP6DE3 open reading frame (ORF) was amplified by PCR using CYP6DE3F1 and CYP6DE3R1 primers (Additional file 1) and CloneAmp HiFi PCR Premix (Takara Bio USA, Inc., Mountain View, CA), cloned into pENTR4 modified to remove the NcoI site , and transformed into Stellar™ Competent Cells (Takara Bio USA, Inc.). Recombinant plasmid was confirmed by sequencing prior to recombination into BaculoDirect™ C-Term Linear DNA (Invitrogen) by LR Clonase™ II (Invitrogen). The recombinant BaculoDirect clone was transferred into Sf9 cells by transfection using Cellfectin II (Invitrogen) and amplified by successive infections of P1 and P2 viral stocks to a high-titer P3 viral stock. Protein expression was initiated by infecting 50 mL of 106 cells/mL Sf9 cells in Sf-900 II SFM culture media supplemented with 10% (vol/vol) fetal bovine serum (Atlas Biologicals, Fort Collins, CO), 0.3 mM δ-aminolevulinic acid, and 0.1 mM ferric citrate with 50 μL of the P3 viral stock and incubating at 27 °C for 72 h. Recombinant CYP6DE3 and housefly cytochrome P450 reductase (HF-CPR)  were harvested 72 h post infection in a cell lysis buffer (100 mM sodium phosphate, pH 7.6, 20% (vol/vol) glycerol, 1.1 mM EDTA, 200 μM PMSF, and protease inhibitor cocktail (Sigma-Aldrich, St. Louis, MO) and the microsomes were isolated by differential centrifugation. The microsomal fraction was tested for functional CYP6DE3 using a CO-difference spectrum analysis .
3-Carene, R-(+)-limonene, (+)-α-pinene and (Z)-dec-7-enal were obtained from Sigma-Aldrich. cis-3-Nonene was obtained from MP Biomedicals (Santa Ana, CA). Reaction mixtures consisted of 200 μL of the CYP6DE3 microsomal fraction, 40 μL of HF-CPR microsomal fraction, 100 mM sodium phosphate buffer pH 7.6, 1.5 mM NADPH (Sigma-Aldrich) and 21 mM of substrate in a total volume of 602 μL. Control reactions containing only HF-CPR were identical to the experimental reactions except that the reaction buffer was substituted for the CYP6DE3 microsomal fraction. Reactions were incubated in a capped 5 mL vial and rotated lengthwise at 30 °C in a FisherBiotech Hybridization Incubator (Thermo Fisher Scientific, Waltham, MA) for three hours. The reactions were terminated and extracted using pentane:ether (1:1). The extracts were analyzed by GC-MS on a HP-5 ms capillary column (Agilent) using an Agilent (Santa Clara, CA) 7890B gas chromatograph coupled to a 5977A mass spectrum detector. The instrument running parameters were: initial temperature of 40 °C with a one min hold, 5 °C/min to 240 °C and 15 °C/min to 300 °C with a 5 min hold. The MS detector was a single quadrapole with an electron ionization source and a molecular weight scanning range of 40 to 700 atomic mass units and an ionization potential of 70 eV.
qRT-PCR of monoterpene exposed beetles
Beetles were separated by sex and placed in two oz plastic cups with perforated lids and incubated in a humidified dark drawer for 24 h, as described previously  to ensure they were unfed at the beginning of the monoterpene exposures. Small clumps of glass wool were placed in four pyrolized glass jars, two of which contained two mL vials capped with 500 μL of a selected monoterpene and a cotton mesh lid, or two jars with no vials as a control. Eight live males or females were transferred into each of the four jars so that each sex had monoterpene-exposed and control treatments. The jars were incubated in the dark for 24 h followed by placing 2-3 beetles in each of three replicate microcentrifuge tubes for each treatment and flash freezing them in liquid nitrogen. Six different monoterpenes were tested: 3-carene, R-(+)-limonene, myrcene, (+)-α-pinene, (–)-α-pinene, terpinolene. We also tested a monoterpene cocktail containing all six listed monoterpenes. RNA was extracted from the whole beetles using the RNeasy Plant Mini Kit from Qiagen as described above with an on-column DNase treatment. qRT-PCR was conducted using CYP6DE3 primers and normalized to YQE_05788 (rpS3P) as described above. Statistically significant differences between the means of relative expression between males and females for each gene were measured using an unpaired two-sample t-test at p < 0.05.
RNA-Seq quality control and validation
Illumina NextSeq500 read processing and mapping results from RNA-seq mountain pine beetle libraries
Total read pairs
Processed read pairs
Adult fed female midgut/fat body
Adult fed female midgut/fat body
Adult fed female midgut/fat body
Adult fed female midgut/fat body
Adult fed male midgut/fat body
Adult fed male midgut/fat body
Adult fed male midgut/fat body
Adult fed male midgut/fat body
Adult unfed female midgut/fat body
Adult unfed female midgut/fat body
Adult unfed female midgut/fat body
Adult unfed female midgut/fat body
Adult unfed male midgut/fat body
Adult unfed male midgut/fat body
Adult unfed male midgut/fat body
Adult unfed male midgut/fat body
Correlation coefficients of RNA-Seq and qRT-PCR measured expression levels across seven samples and 15 genes
petal generated a co-expression network model based on Spearman Correlation Coefficient and similarity threshold of 0.808. This model includes 10,661 transcripts (94% of genes from the entire dataset) and follows the well-established biological network structure characteristics: small-world and scale-free . From this model, closely connected gene groups based upon the genes of interest were extracted, resulting in thirteen vicinity networks (VNs). Here, a vicinity network (VN) is defined by the genes of interests and their common neighbors. For more detail refer to .
Gene ontology enrichment analysis of differential expressed genes
Summary of genes differentially expressed with greater than two fold change and adjusted p-value ≤ 0.01
# of genes
# of genes
Male Fed > Male Unfed
Male Unfed > Male Fed
Male Fed > Female Fed
Female Fed > Male Fed
Male Unfed > Female Unfed
Female Unfed > Male Unfed
Female Fed > Female Unfed
Female Unfed > Female Fed
Pheromone biosynthesis candidate genes
Putative trans-verbenol-biosynthetic genes
Putative exo-brevicomin-biosynthetic genes
exo-Brevicomin production from long chain fatty acid precursors in the fat body of unfed males  involves steps catalyzed by a desaturase, P450s, a short chain dehydrogenase, and a cyclase (Fig. 1c). The mRNA levels of two genes previously identified as active in the biosynthetic pathway, CYP6CR1 and ZnoDH , were elevated in unfed males (Fig. 6b). Other P450 genes with high expression levels in unfed males included CYP6DE3 (YQE_02812), CYP6BW4 (YQE_01441), CYP6DF1 (YQE_11788), CYP4EX1 (YQE_01611), CYP4CV2 (YQE_05823), a short chain dehydrogenase (YQE_04359) and one putative terpenoid cyclase (YQE_04789) (Fig. 6b). The expression profile of YQE_03851, which encodes a P450 with 98% a.a. identity to CYP4G56, similarly had increased expression levels in unfed males.
The purple gene module in Fig. 3a includes seven genes hypothesized or confirmed to be involved in exo-brevicomin biosynthesis (CYP6CR1, ZnoDH, CYP6DE3, CYP4EX1, a cyclase, an oxidoreductase and a hydrolase) and 22 common neighbors. This module is densely connected, with density equal to 0.987 and all genes having greater expression values in males than in females (Fig. 3c, Additional file 3). All gene members are statistically-significantly differentially-expressed (log2-fold change > 1) at a significance level of 0.01 (see Additional file 2).
Putative frontalin-biosynthetic genes
Steps for frontalin biosynthesis downstream of geranylgeranyl diphosphate are likely catalyzed by a dioxygenase or P450, and at least one additional P450 and a cyclase (Fig. 1b). Genes with higher relative expression in fed males were considered as candidates for frontalin biosynthesis. These included four P450s CYP6DK1 (YQE_01078), CYP6DE4 (YQE_01868), CYP345F1 (YQE_06277) and CYP6BW3 (YQE_02884) (Fig. 6c). No putative dioxygenase-encoding genes with elevated transcript levels in fed males compared to other treatments were identified. One putative short chain dehydrogenase, YQE_11963, and one putative terpenoid cyclase, YQE_04789, also had increased mRNA levels in fed males relative to the other groups (Fig. 6b and c).
The light blue gene module in Fig. 3b is based on seven genes previously hypothesized to play a role in frontalin biosynthesis: CYP6DE4, CYP6BW3, a kinase, two oxidoreductases, HMG-CoA reductase (HMGR), and geranylgeranyl diphosphate synthase (GGPPS). These seven genes and their 31 common neighbors created the densely connected light blue module (density = 0.983). The majority (26 of 38 genes) were statistically significantly over-expressed in fed males compared to unfed males and fed females (Fig. 3c).
Another interesting grouping is presented in the orange module. This module also includes seven putative candidate genes for frontalin biosynthesis: CYP345F1, CYP6DK1, CYP6DH2, a CYP6DK1-like P450, a prenyltransferase, an acyl-CoA oxidase, and a dehydrogenase (Fig. 3b). These seven candidate genes and their 77 common immediate neighbors are almost perfectly intra-connected (density = 0.994). Similar to the light blue module, the majority of genes (45) have a statistically significant higher expression (log2-fold change > 1) in fed male than unfed male and fed female (Fig. 3c).
Functional analysis of CYP6DE3
Comparative transcriptome analysis to tentatively identify MPB genes encoding enzymes active in pheromone-biosynthetic pathways is based on the hypothesis that the genes for each pathway are coordinately regulated. It has proven useful for prior studies [15, 18]. We extended this approach here by comparing expression profiles in pheromone-biosynthetic tissues of fed and unfed female and male mountain pine beetles. Extensive RNA-Seq profiling yielded nearly 425 million paired-end reads, with 74.9% aligning to the reference genome. The clear separation of gender and feeding status shown by PCA underscores the remarkable shifts in genome usage exhibited by these beetles [7, 8]. The close correlation between RNA-Seq and qRT-PCR data (Table 2) supports that the expression values reported here reliably indicate in vivo mRNA levels.
We used a combination of bioinformatics analyses to narrow the pool of candidate pheromone-biosynthetic genes, beginning with four straight-forward comparisons of relative expression levels between physiological conditions. For example, exo-brevicomin-biosynthetic genes would be expected to have elevated expression levels in unfed males compared to both unfed females and fed males, whereas frontalin-biosynthetic genes would be higher in fed males compared to both unfed males and fed females. This produced preliminary pools of only ~200 – 1,200 candidates, depending on the analysis (Table 3; Fig. 3). These pools were enriched for P450 production and activity, consistent with increased metabolic activity upon feeding [7, 8, 18], which complicates identifying enzymes involved primarily in pheromone biosynthesis. Nevertheless, mevalonate pathway enzymes are more predominantly represented in fed males compared to females, consistent with frontalin production (Fig. 5). In parallel, co-expression network analysis by petal using a Spearman Correlation Coefficient and similarity threshold of 0.808 also isolated candidate genes. As these expression data were not normally distributed, the Spearman Correlation Coefficient supplied a robust non-parametric alternative to the standard Pearson Correlation Coefficient. Twenty-two of the final genes selected by the petal analysis appear relevant to exo-brevicomin biosynthesis, while another 71 may be involved in frontalin production.
(–)-trans-Verbenol is produced by a single P450-mediated hydroxylation of (–)-α-pinene, a reaction that may be catalyzed by multiple enzymes as part of a detoxification process . Thus, a “pheromone-biosynthetic” P450 that specifically produces trans-verbenol in females may be an artificial designation. Our current study notes three P450s with relatively high expression levels in fed females compared to unfed females and males (Fig. 6a). A fourth P450, CYP4BD4, showed the highest mRNA levels in females compared to males, though in a pattern that is not consistent with feeding-induced trans-verbenol production.
exo-Brevicomin production from long chain fatty acid precursors in the fat body of unfed males  involves steps catalyzed by P450s, a short chain dehydrogenase, and a cyclase (Fig. 1c). High probability candidate genes for exo-brevicomin biosynthesis are likely in the same petal group containing CYP6CR1 and ZnoDH (Fig. 3a). Interestingly, this gene group includes a putative cyclase (YQE_04789) that may catalyze the terminal reaction. The two P450s (CYP6DE3 and CYP4EX1) may be active upstream of ZnoDH to produce and/or hydroxylate 3-nonene. In this respect, the CYP4G56-like P450 (YQE_03851) was not part of the gene group but is of interest given its similar expression profile (Fig. 6b) and identity as a CYP4G. While predicting P450 function from sequences is very difficult, CYP4G family P450s appear to be insect-specific and function as oxidative decarbonylases – yielding hydrocarbons from long chain fatty aldehydes [38, 39]. Thus, YQE_03851 may contribute to 3-nonene production.
Frontalin-biosynthetic steps through the mevalonate pathway to geranylgeranyl diphosphate are well established in fed and JH treated MPB males . Our analysis also identified mRNAs for mevalonate pathway enzymes, including HMGR and GGPPS, to be elevated in fed males compared to other treatment groups. Later steps are likely catalyzed by P450s, a dioxygenase, and a cyclase that should group together with HMGR and GGPPS in the petal analysis. Two P450 genes, CYP6DE4 (YQE_01868) and CYP6BW3 (YQE_02884), did group with HMGR and GGPS (light blue VN in Fig. 3b) while four other P450 genes, CYP345F1 (YQE_06277), CYP6DK1 (YQE_01078), CYP6DH2 (YQE_01329) and a CYP6DK1-like P450 (YQE_01079), grouped into one different VN (orange VN in Fig. 3b). The two VNs are connected directly by two links, and both gene groups portray increased expression in fed males (Fig. 3c), a pattern consistent with frontalin biosynthesis. Interestingly, a putative dioxygenase was not identified, which may suggest alternative activities on a GGPP precursor, perhaps catalyzed by a cytochrome P450. It is also noteworthy that the cyclase identified in the “exo-brevicomin cluster” (YQE_04789) also shows elevated mRNA levels in fed males (Fig. 6b). Given the structural similarities of the epoxide precursors for both exo-brevicomin and frontalin, it is possible that a single cyclase could serve the terminal steps in both pathways.
While comparative transcriptomics is invaluable to preliminarily identify putative pheromone-biosynthetic genes, a more accurate assessment requires additional information . For MPB, our transcriptomic analyses return more candidate genes than there are reactions to catalyze. We hypothesized that those with elevated expression upon exposure to monoterpenes are more likely to contribute to resin detoxification than pheromone component production (except for the case of trans-verbenol, as noted above). We therefore measured relative mRNA levels for CYP6DE3, which we had tagged as a potential exo-brevicomin biosynthetic enzyme, in beetles that had been exposed to atmospheres saturated with various monoterpenes. The clear elevation observed for all cases (Fig. 7) suggests that CYP6DE3 is induced by monoterpene exposure, particularly in females, implying a resin-detoxifying role. The absence of this induction in fed insects further implicates that CYP6DE3 regulation is complex. The monoterpene-dependent difference in response in males and females is curious, but has been exhibited in another study reporting similar sex-specific transcriptional responses of various D. armandi P450 genes in response to monoterpenes . A detoxification role for CYP6DE3 is supported by functional assays of the recombinant enzyme which showed that it oxidized a variety of monoterpenes, but did not appear to accept exo-brevicomin precursors as substrates (Fig. 8 and data not shown). Interestingly, the products at 15.65 min for (+)-α-pinene and 14.78 min for 3-carene have a m/z peak at 168 suggesting these substrates were oxidized twice (Fig. 8a and b).
Candidate genes for MPB pheromone biosynthesis identified by RNA-seq using feeding status or JH treatment
Feeding status and sex
JH treatment and sex (Keeling et al. )
Proposed in both studies
Pheromone biosynthetic pathway
This study identified a number of candidate genes for involvement in MPB aggregation and anti-aggregation pheromone biosynthesis through differential gene expression analysis based on sex and feeding status. However, as Keeling et al.  noted, caution should be employed when using comparative transcriptomic data to identify putative pheromone-biosynthetic genes. As evidenced by the functional analysis of an identified candidate gene, CYP6DE3, for exo-brevicomin biosynthesis in MPB, expression profiles are not always indicative of a specific role for a gene of interest. Further functional analysis of the genes identified in this study should lead to the discovery of most, if not all, of the unknown enzymes involved in MPB aggregation and anti-aggregation pheromone biosynthesis.
Complementary deoxyribonucleic acid
Gas chromatography-mass spectrometry
Georgia Genomics Facility
Geranygeranyl diphosphate synthase
Gene transfer format
House fly cytochrome P450 reductase
Mountain pine beetle
Messenger ribonucleic acid
Nicotinamide adenine dinucleotide phosphate
Open reading frame
Cytochrome P450 monooxygenase
Principal component analysis
Polymerase chain reaction
Quantitative reverse transcriptase PCR
We thank the Georgia Genomics Facility for the Illumina sequencing services and D. Kosma for access to the GC-MS instrument.
The project described was supported by awards from USDA-NIFA (2014-67013-21748), The State of Nevada (NV00381) and National Institute of General Medical Sciences (P20GM103440) from the National Institutes of Health.
Availability of data and materials
The dataset supporting the conclusions of this article is available in the NCBI Sequence Read Archive, https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP072778, at SRA Project Accession SRP072778, with additional project and sample descriptions accessible at NCBI BioProject PRJNA317010, and BioSamples SAMN04595104-SAMN04595119, respectively.
CT, JN, KS and GJB conceived and designed the project. JN, KJ, MF, MM, and SY participated in data acquisition. JN, JP, RJT, and KS analyzed the data. JN, KS, GJB, JP, and CT prepared figures and wrote and edited the manuscript. All authors reviewed and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Mountain Pine Beetles were collected in the Tahoe Donner housing subdivision in Truckee, CA with permission from the Tahoe Donner Forestry Department. No other permissions were necessary.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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.
- Pitman GB, Vité JP, Kinzer GW, Fentiman AF. Bark Beetle Attractants: Trans-verbenol isolated from Dendroctonus. Nature. 1968;218:168–9.View ArticleGoogle Scholar
- Song M, Delaplain P, Nguyen TT, Liu X, Wickenberg L, Jeffrey C, et al. exo-Brevicomin biosynthetic pathway enzymes from the Mountain Pine Beetle, Dendroctonus ponderosae. Insect Biochem Mol Biol. 2014;53:73–80.View ArticlePubMedGoogle Scholar
- Song M, Gorzalski A, Nguyen TT, Liu X, Jeffrey C, Blomquist GJ, et al. exo-Brevicomin biosynthesis in the fat body of the mountain pine beetle, Dendroctonus ponderosae. J Chem Ecol. 2014;40:181–9.View ArticlePubMedGoogle Scholar
- Pureswaran DS, Gries R, Borden JH, Pierce Jr HD. Dynamics of pheromone production and communication in the mountain pine beetle, Dendroctonus ponderosae Hopkins, and the pine engraver, Ips pini (Say) (Coleoptera: Scolytidae). Chemoecology. 2000;10:153–68.View ArticleGoogle Scholar
- Barkawi LS, Francke W, Blomquist GJ, Seybold SJ. Frontalin: De novo biosynthesis of an aggregation pheromone component by Dendroctonus spp. bark beetles (Coleoptera: Scolytidae). Insect Biochem Mol Biol. 2003;33:773–88.View ArticlePubMedGoogle Scholar
- Keeling CI, Yuen MM, Liao NY, Roderick Docking T, Chan SK, Taylor GA, et al. Draft genome of the mountain pine beetle, Dendroctonus ponderosae Hopkins, a major forest pest. Genome Biol. 2013;14:R27.View ArticlePubMedPubMed CentralGoogle Scholar
- Keeling CI, Henderson H, Li M, Yuen M, Clark EL, Fraser JD, et al. Transcriptome and full-length cDNA resources for the mountain pine beetle, Dendroctonus ponderosae Hopkins, a major insect pest of pine forests. Insect Biochem Mol Biol. 2012;42:525–36.View ArticlePubMedGoogle Scholar
- Robert JA, Pitt C, Bonnett TR, Yuen MMS, Keeling CI, Bohlmann J, et al. Disentangling detoxification: gene expression analysis of feeding mountain pine beetle illuminates molecular-level host chemical defense detoxification mechanisms. PLoS One. 2013;8:e77777. Public Library of Science.View ArticlePubMedPubMed CentralGoogle Scholar
- Pitt C, Robert JA, Bonnett TR, Keeling CI, Bohlmann J, Huber DPW. Proteomics indicators of the rapidly shifting physiology from whole mountain pine beetle, Dendroctonus ponderosae (Coleoptera: Curculionidae), adults during early host colonization. PLoS One. 2014;9, e110673.View ArticlePubMedPubMed CentralGoogle Scholar
- Blomquist GJ, Figueroa-Teran R, Aw M, Song M, Gorzalski A, Abbott NL, et al. Pheromone production in bark beetles. Insect Biochem Mol Biol. 2010;40:699–712.View ArticlePubMedGoogle Scholar
- Keeling CI, Li M, Sandhu HK, Henderson H, Yuen MM S, Bohlmann J. Quantitative metabolome, proteome and transcriptome analysis of midgut and fat body tissues in the mountain pine beetle, Dendroctonus ponderosae Hopkins, and insights into pheromone biosynthesis. Insect Biochem Mol Biol. 2016;70:170–83.View ArticlePubMedGoogle Scholar
- Borden JH, Nair KK, Slater CE. Synthetic Juvenile Hormone: Induction of Sex Pheromone Production in Ips confusus. Science. 1969;166:1626–7.View ArticleGoogle Scholar
- Chen NM, Borden JH, Pierce HD. Effect of juvenile hormone analog, fenoxycarb, on pheromone production byIps paraconfusus (Coleoptera: Scolytidae). J Chem Ecol. 1988;14:1087–98.View ArticlePubMedGoogle Scholar
- Tillman JA, Holbrook GL, Dallara PL, Schal C, Wood DL, Blomquist GJ, et al. Endocrine regulation of de novo aggregation pheromone biosynthesis in the pine engraver, Ips pini (Say) (Coleoptera: Scolytidae). Insect Biochem Mol Biol. 1998;28:705–15. United Kingdom.View ArticleGoogle Scholar
- Keeling CI, Bearfield JC, Young S, Blomquist GJ, Tittiger C. Effects of juvenile hormone on gene expression in the pheromone-producing midgut of the pine engraver beetle, Ips pini. Insect Mol Biol. 2006;15:207–16.View ArticlePubMedGoogle Scholar
- Tittiger C, Keeling CI, Blomquist GJ. Some Insights into the Remarkable Metabolism of the Bark Beetle Midgut. Recent Adv Phytochem Elsevier. 2005;39:57–78.View ArticleGoogle Scholar
- Bearfield JC, Henry AG, Tittiger C, Blomquist GJ, Ginzel MD. Two regulatory mechanisms of monoterpenoid pheromone production in Ips spp. of bark beetles. J Chem Ecol. 2009;35:689–97.View ArticlePubMedGoogle Scholar
- Aw T, Schlauch K, Keeling CI, Young S, Bearfield JC, Blomquist GJ, et al. Functional genomics of mountain pine beetle (Dendroctonus ponderosae) midguts and fat bodies. BMC Genomics. 2010;11:215. BioMed Central.View ArticlePubMedPubMed CentralGoogle Scholar
- Fleige S, Pfaffl MW. RNA integrity and the effect on the real-time qRT-PCR performance. Mol Aspects Med. 2006;27:126–39.View ArticlePubMedGoogle Scholar
- Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet journal. 2011;17:10.View ArticleGoogle Scholar
- Zimin AV, Marçais G, Puiu D, Roberts M, Salzberg SL, Yorke JA. The MaSuRCA genome assembler. Bioinformatics. 2013;29:2669–77.View ArticlePubMedPubMed CentralGoogle Scholar
- Kersey PJ, Allen JE, Armean I, Boddu S, Bolt BJ, Carvalho-Silva D, et al. Ensembl Genomes 2016: more genomes, more complexity. Nucleic Acids Res. 2015;44:D574–80.View ArticlePubMedPubMed CentralGoogle Scholar
- Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.View ArticlePubMedPubMed CentralGoogle Scholar
- Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.View ArticlePubMedGoogle Scholar
- Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.View ArticlePubMedPubMed CentralGoogle Scholar
- Mitchell A, Chang H-Y, Daugherty L, Fraser M, Hunter S, Lopez R, et al. The InterPro protein families database: the classification resource after 15 years. Nucleic Acids Res. 2014;43:D213–21.View ArticlePubMedPubMed CentralGoogle Scholar
- Kinsella RJ, Kähäri A, Haider S, Zamora J, Proctor G, Spudich G, et al. Ensembl BioMarts: a hub for data retrieval across taxonomic space. Database (Oxford). 2011;2011:bar030.View ArticleGoogle Scholar
- Toutges MJ, Hartzer K, Lord J, Oppert B. Evaluation of reference genes for quantitative polymerase chain reaction across life cycle stages and tissue types of Tribolium castaneum. J Agric Food Chem. 2010;58:8948–51.View ArticlePubMedGoogle Scholar
- Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25:402–8.View ArticlePubMedGoogle Scholar
- Dehmer M, Emmert-Streib F, Graber A, Salvador A. Applied Statistics for Network Biology: Methods in Systems Biology. Germany: Wiley-Blackwell; 2011.View ArticleGoogle Scholar
- Petereit J, Smith S, Harris FC, Schlauch KA. petal: Co-expression network modelling in R. BMC Syst Biol. 2016;10 Suppl 2:51.View ArticlePubMedPubMed CentralGoogle Scholar
- Benjamini Y, Hochberg Y. Controlling The False Discovery Rate - A Practical And Powerful Approach To Multiple Testing. J R Stat Soc Ser B Methodol. 1995;57:289–300.Google Scholar
- Maere S, Heymans K, Kuiper M. BiNGO: a Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics. 2005;21:3448–9.View ArticlePubMedGoogle Scholar
- Sandstrom P, Welch WH, Blomquist GJ, Tittiger C. Functional expression of a bark beetle cytochrome P450 that hydroxylates myrcene to ipsdienol. Insect Biochem Mol Biol. 2006;36:835–45.View ArticlePubMedGoogle Scholar
- Choi S-J, Kim M, Kim S-I, Jeon J-K. Microplate assay measurement of cytochrome p450-carbon monoxide complexes. J Biochem Mol Biol. 2003;36:332–5.PubMedGoogle Scholar
- Pierce HD, Conn JE, Oehlschlager AC, Borden JH. Monoterpene metabolism in female mountain pine beetles, Dendroctonus ponderosae Hopkins, attacking ponderosa pine. J Chem Ecol. 1987;13:1455–80.View ArticlePubMedGoogle Scholar
- Qiu Y, Tittiger C, Wicker-Thomas C, Le Goff G, Young S, Wajnberg E, et al. An insect-specific P450 oxidative decarbonylase for cuticular hydrocarbon biosynthesis. Proc Natl Acad Sci U S A. 2012;109:14858–63.View ArticlePubMedPubMed CentralGoogle Scholar
- Balabanidou V, Kampouraki A, MacLean M, Blomquist GJ, Tittiger C, Juárez MP, et al. Cytochrome P450 associated with insecticide resistance catalyzes cuticular hydrocarbon production in Anopheles gambiae. Proc Natl Acad Sci U S A. 2016;113:9268–73.View ArticlePubMedPubMed CentralGoogle Scholar
- Keeling CI, Chiu CC, Aw T, Li M, Henderson H, Tittiger C, et al. Frontalin pheromone biosynthesis in the mountain pine beetle, Dendroctonus ponderosae, and the role of isoprenyl diphosphate synthases. Proc Natl Acad Sci U S A. 2013;110:18838–43.View ArticlePubMedPubMed CentralGoogle Scholar
- Dai L, Ma M, Wang C, Shi Q, Zhang R, Chen H. Cytochrome P450s from the Chinese white pine beetle, Dendroctonus armandi (Curculionidae: Scolytinae): Expression profiles of different stages and responses to host allelochemicals. Insect Biochem Mol Biol. 2015;65:35–46.View ArticlePubMedGoogle Scholar
- Mock KE, Bentz BJ, O’neill EM, Chong JP, Orwin J, Pfrender ME. Landscape-scale genetic variation in a forest outbreak species, the mountain pine beetle (Dendroctonus ponderosae). Mol Ecol. 2007;16:553–68.View ArticlePubMedGoogle Scholar