Comparative transcriptomics of mountain pine beetle pheromone-biosynthetic tissues and functional analysis of CYP6DE3
BMC Genomics volume 18, Article number: 311 (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.
The mountain pine beetle (MPB) (Dendroctonus ponderosae Hopkins) uses three main pheromone components to coordinate the “mass attack” necessary to overcome a host tree’s defenses. Each component has a distinctive role and is produced from a different metabolic pathway (Fig. 1). Females produce the aggregation pheromone (–)-trans-verbenol [(1S,2R,5S)-4,6,6-trimethylbicyclo[3.1.1] hept-3-en-2-ol], to attract other males and females to a newly colonized host tree . Biosynthesis of (–)-trans-verbenol likely requires cytochrome P450-mediated hydroxylation of the host tree produced monoterpene, (–)-α-pinene. Newly emerged males produce the aggregation pheromone exo-brevicomin [(1R,5S,7R)-7-ethyl-5-methyl-6,8-dioxabicyclo[3.2.1]octane] from long chain fatty acid precursors in the fat body [2, 3]. exo-Brevicomin production decreases substantially upon reaching the new host tree and mating . The third major component, frontalin [(1R,5R)-1,5-dimethyl-6,8-dioxabicyclo[3.2.1]octane], is produced by feeding males and acts as an anti-aggregation signal to halt the attack and prevent overcrowding of the tree . Frontalin biosynthesis uses the mevalonate pathway .
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
A total of 424,776,657 paired-end reads consisting of at least 76 bp were recovered from the 16 libraries, with reads per library ranging from 18,659,429 to 33,943,439 with a mean of 26,548,541. After processing, the number of reads aligned to the reference genome was 317,944,928 (Table 1), representing 13,407 potentially different transcripts. Of these transcripts a total of 11,342 transcripts passed the two filtering steps as explained in the Materials and Methods. Verification of RNA-Seq by qRT-PCR analysis confirmed similar expression measures between the two platforms. For the 15 genes in the seven samples measured by qRT-PCR and RNA-Seq, Pearson and Spearman Correlation Coefficients averaged to 0.924 and 0.878, respectively (Table 2). The PCA showed a clear separation between feeding states and genders, with almost 70% of the variance explained by feeding states (Fig. 2).
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 .
Based on the results from petal, further in-depth analyses were conducted leading to three distinct gene modules, one representing candidates involved in exo-brevicomin biosynthesis (purple), and two representing candidates involved in frontalin biosynthesis (light blue and orange) (Fig. 3a and b). Gene membership of each module is listed in Additional file 3.
Gene ontology enrichment analysis of differential expressed genes
Based on the differential gene analysis 3,799 differentially expressed transcripts had a statistically significant (log2-fold change > 1) differential gene expression at a significance level of 0.01 after multiple testing adjustment (Table 3, Additional file 4). The greatest number of genes with a differential expression matching these criteria occurred in the MU > FU, MU > MF, and MF > MU comparisons with 894, 698, and 638 genes, respectively. The smallest numbers of genes with a differential expression matching our criteria were in the comparisons of FU > MU and FF > MF with 74 and 108 genes respectively.
Two areas highlighted in the Venn diagram (Fig. 4) were of particular interest as they represented pools of potential candidate genes for pheromone biosynthesis: 1) the 1225 genes with statistically significant differential expression between MU vs. FU beetles are exo-brevicomin-biosynthetic candidates; 2) the 217 genes with statistically significant differential expression in MF vs. MU and MF vs. FF, excluding differentially expressed genes in MU vs. FU, are potential frontalin-biosynthetic candidates.
To confirm potential candidate genes for pheromone biosynthesis, GO enrichment analysis was performed. Of the final 11,342 transcripts, 7,603 had at least one GO term. For the two pools of statistically significant differentially-expressed genes identified in the Venn diagram (possible exo-brevicomin- and frontalin-biosynthetic genes), a GO enrichment analysis using BiNGO within Cytoscape identified over-represented molecular functions and metabolic processes with an adjusted significance value of less than 0.05 (Fig. 5). Enriched GO terms in the possible exo-brevicomin-biosynthetic genes included P450-associated terms such as iron ion binding, heme binding and increased monooxygenase activity. Within the possible frontalin-biosynthetic genes, a wider variety of enriched metabolic processes are represented: uroporphyrinogen-III synthase and tetrapyrrole biosynthetic activity, monoxygenase activity, carbohydrate metabolic processes, isoprenoid biosynthetic processes, and two processes involved in P450 biosynthesis. Both gene groups showed enrichment in P450 related processes, however GO terms related to the mevalonate pathway were more highly enriched in the putative frontalin-synthesizing candidates.
Pheromone biosynthesis candidate genes
Putative trans-verbenol-biosynthetic genes
trans-Verbenol biosynthesis in female MPB likely takes place via P450-mediated hydroxylation of host (–)-α-pinene . Three out of four P450 genes, CYP6DJ1 (YQE_04799), CYP6DJ2 (YQE_04800) and CYP349B2 (YQE_02158) showed relatively high mRNA levels in fed females (Fig. 6a). One other P450, CYP4BD4 (YQE_07200), showed mRNA levels in both fed and unfed females higher than those in males (Fig. 6a).
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
CYP6DE3 mRNA levels were elevated in all samples of unfed male and female beetles exposed to a variety of monoterpenes for 24 h compared to the controls. In general, monoterpene exposure elevated CYP6DE3 mRNA levels in females more strongly than in males, with the exception that (+)-α-pinene exposure elevated mRNA levels more strongly in males (Fig. 7).
Recombinant CYP6DE3 was assayed for activity with hypothesized intermediates in the exo-brevicomin biosynthetic pathway, (Z)-dec-7-enal and (Z)-3-nonene (Fig. 1), as well as a variety of monoterpenes. Assays with (Z)-dec-7-enal or (Z)-3-nonene showed no unique products when analyzed by GC-MS compared to the negative controls, and neither were (Z)-3-nonene nor 6(Z)-nonen-2-ol, the hypothesized products, respectively, detected. The monoterpene substrates, however, showed unique or substantially increased GC-MS retention peaks compared to controls. Products resulting from (+)-α-pinene were tentatively identified as 3-oxatricyclo [220.127.116.11(2, 4)] octane (10.88 min) and trans-verbenol (12.15 min). Two other products (at 15.42 and 15.65 min) remain unidentified (Fig. 8a). Similarly, incubations with 3-carene as a substrate yielded two unknown products at 14.78 and 15.92 min (Fig. 8b), and incubations with (+)-limonene as a substrate also produced two unknown products (16.00 and 16.12 min; Fig. 8c).
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).
De novo pheromone component biosynthesis in pine bark beetles is affected by sex, feeding status, environment, and JH III , with JH III treatment sometimes being sufficient to elevate mRNA levels of pheromone-biosynthetic genes even in insects that otherwise require feeding to trigger pheromone production . Indeed, JH III stimulates both frontalin [5, 11] and trans-verbenol biosynthesis, but not exo-brevicomin biosynthesis  in MPB. Our study complements those of Robert et al. , who compared fed and JH III-treated whole insects and concentrated on a survey of detoxification mechanisms, and Keeling et al. , who compared starved and JH III-treated midguts and fat bodies. Our study differs in that we focused on midgut and fat body tissues of fed and unfed insects rather than JH III-treated insects because of the evident complexity in regulating production of these three main pheromone components. Several putative pheromone-biosynthetic genes identified in our study agree with those reported by Keeling et al.  (Table 4), and the increased confidence resulting from this concurrence makes the common genes high priorities for functional assays. It is also noteworthy that CYP6DE4 does not accept pheromone precursors despite being induced by JH III . The discrepancies in the list of candidate enzymes are likely due to a combination of factors, including differences in experimental design and data analysis. Given that the populations used by Keeling et al.  and us appear to be geographically and genetically isolated , it is also possible that their responses to different conditions also differ (, unpublished data).
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
Pitman GB, Vité JP, Kinzer GW, Fentiman AF. Bark Beetle Attractants: Trans-verbenol isolated from Dendroctonus. Nature. 1968;218:168–9.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Borden JH, Nair KK, Slater CE. Synthetic Juvenile Hormone: Induction of Sex Pheromone Production in Ips confusus. Science. 1969;166:1626–7.
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.
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.
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.
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.
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.
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.
Fleige S, Pfaffl MW. RNA integrity and the effect on the real-time qRT-PCR performance. Mol Aspects Med. 2006;27:126–39.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet journal. 2011;17:10.
Zimin AV, Marçais G, Puiu D, Roberts M, Salzberg SL, Yorke JA. The MaSuRCA genome assembler. Bioinformatics. 2013;29:2669–77.
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.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.
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.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
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.
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.
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.
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.
Dehmer M, Emmert-Streib F, Graber A, Salvador A. Applied Statistics for Network Biology: Methods in Systems Biology. Germany: Wiley-Blackwell; 2011.
Petereit J, Smith S, Harris FC, Schlauch KA. petal: Co-expression network modelling in R. BMC Syst Biol. 2016;10 Suppl 2:51.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Primer Table: A list of primers used for qRT-PCR and CYP6DE3 cloning. (XLSX 10 kb)
petal Venn Diagrams: Venn diagrams showing the intersections of differentially over-expressed genes with a < 0.01 for the (A) purple, (B) light blue and (C) orange gene groups identified in the petal analysis. “gr” indicates greater than in each of the comparisons. (PPTX 178 kb)
petal Module Gene Annotations: A list of genes and their annotations found in each petal module. (XLSX 20 kb)
Differentially Expressed Genes for Each Comparison: A list of all the genes that have greater than two fold differential expression for each considered comparison of sex and feeding status. (XLSX 1704 kb)
About this article
Cite this article
Nadeau, J.A., Petereit, J., Tillett, R.L. et al. Comparative transcriptomics of mountain pine beetle pheromone-biosynthetic tissues and functional analysis of CYP6DE3. BMC Genomics 18, 311 (2017). https://doi.org/10.1186/s12864-017-3696-4