Introduction
============

*Pseudomonas* spp. and their genomic features (e.g., multifaceted oxidative enzymatic repertoire) have been extensively explored as cell-factories and functional blocks for biocatalysis ([@evu254-B85]; [@evu254-B83]; [@evu254-B18]; [@evu254-B69]; [@evu254-B59]), particularly in the production of fine chemicals for several industries (e.g., terpenoids) ([@evu254-B13]; [@evu254-B55]; [@evu254-B79]; [@evu254-B21]; [@evu254-B23]; [@evu254-B9]). The instability of some monoterpenes directed biotechnological research to their derivatives which are known to retain or enhance some properties of the original compound ([@evu254-B79]) and are usually easier to incorporate into hydrophilic/aqueous formulations.

The utilization of monoterpenes and their derivatives is widespread in industry because of their unique properties ([@evu254-B37]; [@evu254-B17]; [@evu254-B61]). One of the most promising monoterpenes, displaying a plethora of industrial applications, is β-myrcene (7-Methyl-3-methylene-1,6-octadiene). β-myrcene is mainly used for the synthesis of top selling flavors and fragrances (e.g., linalool, nerol, geraniol, menthol) ([@evu254-B8]) and compounds with pharmacological potential (e.g., antimutagenics, analgesics, tyrosinase inhibitors), as well as a starting material in production of polymers, biodegradable surfactants, pheromones, and insect repellents ([@evu254-B43]; [@evu254-B53]; [@evu254-B8]).

In 1999, a *Pseudomonas* sp. (strain M1) was isolated from soil sediments of the Rhine River ([@evu254-B40]), which exhibited the ability to use β-myrcene as sole carbon and energy sources. Transposon random mutagenesis led to the identification of four genes potentially coding for β-myrcene-biotransforming enzymes: An aldehyde dehydrogenase (*myrA*), an alcohol dehydrogenase (*myrB*), an acyl-CoA synthetase (*myrC*), and an enoyl-CoA hydratase (*myrD*) ([@evu254-B40]). Later, a proteomic analysis ([@evu254-B68]) provided further insights into the β-myrcene catabolic pathway and associated adaptive mechanisms of *Pseudomonas* sp. M1. However, during that proteomic analysis, the authors outlined that a significant number of β-myrcene-responsive proteins were not identified due to the lack of information in the public databases, suggesting that the β-myrcene catabolic apparatus of *Pseudomonas* sp. M1 would be based, at least partially, on enzymes that have evolved divergently with respect to what have been reported in the literature so far. Curiously, in another work using *Pseudomonas* sp. M1, it was reported that its ability to use phenol and benzene as sole carbon sources relied on a unique σ^54^-dependent gene cluster organization ([@evu254-B67]), which, apparently, has evolved differently, when compared with other reported phenol catabolic pathways. Altogether, these hints suggest that *Pseudomonas* sp. M1 harbors, in its genomic repertoire, a unique code for biotechnologically relevant enzymes and functional modules, probably due to particular evolutionary constrains.

Nowadays, *Pseudomonas* sp. M1 is one of the most promising gram-negative bacteria able to metabolize β-myrcene, for which a catabolic pathway draft has been proposed ([@evu254-B40]; [@evu254-B68]). *Pseudomonas* sp. M1 is also capable of mineralizing other terpenes (e.g., citronellol, citral), as well as several toxic and/or recalcitrant phenolic compounds, such as phenol, 4-chlorophenol, benzene, and toluene ([@evu254-B66]; [@evu254-B68]).

The genome reconstruction conducted in this work updated the current available draft reported by [@evu254-B72] for exploitation of *Pseudomonas* sp. M1 full potential. A genomic analysis led to the identification of genomic loci that suggest a versatile metabolism of terpene-backbone compounds. Transcriptome sequencing of β-myrcene-challenged M1 cells revealed that the core-code of β-myrcene metabolism is located in a unique genetic locus and allowed the identification of the full set of M1 genes involved in β-myrcene-dependent molecular mechanisms of catabolism, adaption, and regulation.

Materials and Methods
=====================

Bacterial Strains, Culture Conditions, and Sampling
---------------------------------------------------

*Pseudomonas* sp. M1 ([@evu254-B40]) was grown in Lysogenic Broth Lennox for DNA extraction. For gene expression analysis, cells were cultivated in 50 ml of mineral medium (MM; 8.9 mM K~2~HPO~4~, 6.2 mM NaH~2~PO~4~, 34.2 μM EDTA, 7 μM ZnSO~4~, 6.8 μM CaCl~s~, 18 μM FeSO~4~, 0.8 μM Na~s~MoO~4~, 0.7 μM CuSO~4~, 1.7 μM CoCl~2~, 1.9 μM MnCl~2~, 15.1 mM (NH~4~)~2~SO~4~, and 0.5 mM MgCl~2~) supplemented with 0.4% lactate or 100 μl of β-myrcene (Merck) (CAS number: 123-35-3) as sole carbon source, in 250-ml Erlenmeyer-sealed flasks, at 30 °C and with orbital shaking (200 rpm). As the water solubility of β-myrcene is very low (density: 0.796 g/ml; log K~ow~: 4.17; water solubility at 25 °C approximately 6 mg/l ([@evu254-B31]), the amount (100 μl) of β-myrcene used in the experiments assured a constant water phase saturation with β-myrcene to avoid eventual carbon starvation effects being reflected on the overall transcriptome analysis.

For total RNA sample preparation, after overnight growth, M1 cells were inoculated to an initial OD~600nm~ of 0.08--0.1 in MM supplemented with respective carbon sources (0.4% lactate or 100 μl of β-myrcene) and cultures diluted again to OD~600nm~ of 0.08--0.1 when OD~600nm~ of 0.5 was reached to synchronize cell growth. Samples for RNA extraction were collected in triplicate (lactate: L1, L2, L3; β-myrcene: M1, M2, M3) when fresh cultures reached an OD~600nm~ of about 0.3--0.4 (2 h of exponential growth) and were immediately flash-frozen in liquid nitrogen.

Next-Generation Genome Sequencing, Assembly, and Annotation
-----------------------------------------------------------

The genomic DNA from *Pseudomonas* sp. M1 was extracted and purified using the Wizard Genomic purification kit (Promega). The DNA samples were processed according to Illumina instructions generating Nextera XT paired-end libraries (2 × 250 bp) which were sequenced using the Illumina MiSeq system at Yale Center for Genomic Analysis.

Raw paired-reads obtained from Illumina MiSeq platform (2 × 250 bp with insert size around 470 bp) and the previously obtained data set using Illumina GAIIx (2 × 50 bp paired-end reads with insert size around 320 bp) ([@evu254-B72]) were both corrected with QUAKE ([@evu254-B44]). These QUAKE-corrected read-data sets and the previously obtained data sets: Illumina GAIIx 2 × 50 bp mate-pair reads with insert size around 5,200 bp, and 454 FLX (single reads with average length of 523 bp) technologies ([@evu254-B72]) were then all processed using fastq-mcf tool ([@evu254-B5]), with which a quality filter was applied on all reads based on Phred quality scores ([@evu254-B20]) and eventual adapter contamination, low quality, and ambiguous nucleotides were trimmed off from the remaining reads. A comparison of the overall quality of the resulting corrected and filtered data sets with the original data sets was carried out using FastQC ([@evu254-B2]).

Paired-end reads from MiSeq were de novo assembled using Edena version 3.131028 ([@evu254-B34]) and Spades version 2.5.1 ([@evu254-B6]) using different parameterization. To analyze the assembly quality, internal SNVs and indels were estimated from reference assembly mapping of the filtered Illumina GAIIx paired-end reads against the respective de novo genome assemblies, using bowtie2 ([@evu254-B45]). Samtools version 0.1.19 ([@evu254-B46]) and SNVer version 0.5.1 ([@evu254-B82]) were used to convert .SAM files to sorted .BAM files and detection of variant nucleotides, respectively. The tuned assemblies from Edena and Spades assemblers were merged using GAM-NGS assembly merger ([@evu254-B80]), followed by automatically improvement of the resulting de novo genome assembly, using Pilon version 1.6 (<http://www.broadinstitute.org/software/pilon/>, last accessed October 1, 2014). Afterwards, filtered data sets from Illumina MiSeq and GAIIx (both) were used for scaffolding with SSPACE version 2.3 ([@evu254-B10]) and gap closure by using Gapfiller version 1.5 ([@evu254-B11]) and GapCloser version 1.12 ([@evu254-B50]).

The final genome sequence draft of *Pseudomonas* sp. M1 was obtained following automatic improvement using Pilon, mapping of Sanger sequencing reads of regions evidencing abnormal coverage and/or gaps, mapping of 454 Roche FLX data set and manual curation. Genome annotation was carried out using National Center for Biotechnology Information (NCBI)'s Prokaryotic Genomes Annotation Pipeline 2.0 ([@evu254-B3]) and metabolic pathways were predicted in silico by KAAS ([@evu254-B54]).

*Total RNA Extraction from Pseudomonas*sp. *M1*
-----------------------------------------------

Total RNA samples were obtained using the Aurum Total RNA Mini kit (Bio-Rad), as described by the manufacturer's instructions. Prior to RNA sequencing (RNA-seq) and reverse transcription-quantitative polymerase chain reaction (RT-qPCR) experiments, all genomic DNA remaining in RNA samples was digested with RNA-free DNase I. RNA integrity was verified in the Agilent 2100 Bioanalyzer (Agilent Technologies, CA) and only samples with RNA Integrity Number over 8.0 were used in RNA-seq and RT-qPCR experiments.

Differential Expression Analysis of RNA-Seq Samples from M1 Cells Grown in Lactate or β-Myrcene
-----------------------------------------------------------------------------------------------

Total RNA samples extracted from M1 cells were processed according to Illumina instructions generating libraries of 50-bp single reads, which were sequenced by Illumina GAIIx platform at BaseClear (the Netherlands), resulting in three data sets from M1 cells grown in lactate medium and three data sets from M1 cells grown in β-myrcene medium, which ranged from 20 million to 56 million reads.

Raw reads were trimmed and clipped for adaptor sequences and ambiguous nucleotides (N's) using fastq-mcf tool ([@evu254-B5]). Expression levels were estimated by mapping and counting filtered RNA-seq data sets against the M1 draft genome using EDGE-pro pipeline ([@evu254-B51]). Differential expression levels were determined as fold change (FC) by the ratio of normalized β-myrcene-derived read count values per lactate values using the EdgeR Exact Test ([@evu254-B63]) with false discovery rate (FDR) *P* value correction. Genes were defined as significantly differentially expressed if they evidenced an altered FC ratio of 1.5 (or higher) and FDR smaller than 0.05. Genes having less than five reads mapped were considered to not be expressed.

*Characterization of the 28-kb Genomic Island Harboring β-Myrcene Core-Code*

In silico prediction of the transcriptional units (TUs), promoter and terminator regions was carried out by FGENESB tool from Softberry server (<http://www.softberry.com/>, last accessed October 1, 2014). RNA-seq based prediction of operons and transcription start sites (TSS) within this genomic island (GI) was carried out using ReadXplorer ([@evu254-B36]).

Quantitative Expression of β-Myrcene-Biotransforming Key Genes by RT-qPCR
-------------------------------------------------------------------------

For RT-qPCR experiments, reverse transcription of RNA to cDNA was performed with the iScript cDNA Synthesis Kit (Bio-Rad). A total of 500 ng of RNA (from the same batch used for RNA-seq) was used for each sample and cDNA synthesized according to manufacturer's instructions (5 min at 25 °C, 30 min at 40 °C, and 5 min at 85 °C).

Expression of *PM1_0216305*, *PM1_0216325*, *PM1_0216335*, *PM1_0216350*, *myrA* (*PM1_0216395*), *myrD* (*PM1_0216400*), and *recA* (*PM1_0217885*) was analyzed in a CFX96 Touch Real-Time PCR Detection System (Bio-Rad) using SensiFAST SYBR No-ROX Kit (Bioline) as the detecting agent. Primers used were synthesized by Stab Vida with following sequences: *PM1_0216305* (fwd: 5′-TATCGCAGCTTGGTCTTGGC-3′; rev: 5′-CTAACCAGTCACCAGACGCC-3′), *PM1_0216325* (fwd: 5′-AAAGGACTGGGCAAACTGGG-3′; rev: 5′-GCTGACGGCTTTGTCTTTGG-3′), *PM1_0216335* (fwd: 5′-CCGGATTACCGAGCTTCTGG-3′; rev: 5′-GGTGTTGGCTGTGTCAGGC-3′), *PM1_0216350* (fwd: 5′-TGTTCTCAAGGTCGCTGACG-3′; rev: 5′-TCCTGCAAATGCTTCAGGGG-3′), *myrA* (fwd: 5′-CGATACCGCTCTGCAACTCC-3′; rev: 5′-CCGCTGCGATGATTTTGTCG-3′), *myrD* (fwd: 5′-GGCGACTTACAGCGTTTTGC-3′; rev: 5′-GCTGGAGTGGTATGCAGAGC-3′), and *recA* (fwd: 5′-GGCCGAGTTCCAGATCCTCTAC-3′; rev: 5′-GATTTCCGCGTTCTCTTCCAGG-3′). Efficiency of each pair of primers (E) was determined with cDNA dilution curves (from 1 μg to 10 pg) from M1 cultures grown in lactate or β-myrcene, using the equation E = 10^(−1/slope)^ − 1 ([@evu254-B65]). Primers used showed efficiencies between 90% and 110%.

Each reaction was set up in a total volume of 20 μl, containing 100 ng of cDNA, SensiFAST SYBR No-ROX mix (Bio-Rad), 2 pmol of each gene-specific primer, and nuclease-free water. Amplification was performed at 95 °C for 2 min followed by 40 cycles at 95 °C for 5 s, 55 °C for 10 s, and 72 °C for 20 s. Specificity of primer annealing and presence of amplification artifacts were assessed by the dissociation temperatures of amplification products in a melting curve analysis of temperatures between 65 °C and 95 °C (stepwise increase of 0.5 °C/5 s), with a continuous measurement of fluorescence. Negative control without cDNA template and control of genomic DNA contamination in RNA samples were performed for each pair of primers and each batch of RNA extracted. Eight independent RNA samples were used and all reactions were performed in triplicated.

Quantification cycle values (Cq) were determined with Bio-Rad CFX Manager software as the threshold cycle. Relative quantification of gene expression as FC was calculated by the 2^−ΔΔCq^ method taking into account primer efficiencies, as described by [@evu254-B58] and using *recA* (*PM1_02298*) as reference gene for normalization ([@evu254-B52]; [@evu254-B29]). Error propagation was calculated according to [@evu254-B57].

Phylogeny of Genes from the 28-kb GI
------------------------------------

Phylogenetic analysis was carried for all 22 genes comprised in the 28-kb GI and for the *GMP synthase* (*guaA*, *PM1_216145*). Each gene was compared with sequences deposited in GenBank by BLASTn ([@evu254-B1]) and the 20 most homologous entries were chosen for multialignment with MAFFT v7 ([@evu254-B42]). Using these alignments, maximum likelihood (ML) trees were inferred with PHYML v3.0 ([@evu254-B30]), using aproximate likelihood-ratio test nonparametric branch support based on a Shimodaira--Hasegawa-like procedure. For all the genes of the GI, the preferred model of nucleotide substitution was GTR (General Time Reversible) with a proportion of invariable sites. Similarities/dissimilarities between phylogenetic tree topologies may provide important clues on the relatedness of the evolutionary track of each gene of the GI under analysis in this work. Therefore, the phylogenetic tree-based parameters Log-likelihood score and patristic distance between multialigned sequences for each gene were deduced from the phylogenetic tree reconstruction and used to build a multivariate matrix. Hierarchical clustering analysis of the resulting matrix was carried by estimating the squared Euclidean distance between the different objects (gene alignment tree parameters), followed by agglomerative clustering with UPGMA (Unweighted Pair Group Method with Arithmetic Mean) method, using R statistical software with the gplots package ([@evu254-B81]). The number of clusters was estimated by partitioning around mediods clustering algorithm, using the function pamk() in the R package fpc ([@evu254-B33]).

Results and Discussion
======================

Whole-Genome Sequencing and Assembly
------------------------------------

Previous work describing key players of β-myrcene catabolism in *Pseudomonas* sp. M1 ([@evu254-B40]; [@evu254-B68]) was limited by the absence of reference genomic data. To set the proper background for exploring the biotransformation potential of *Pseudomonas* sp. M1, its genome was resequenced with Illumina MiSeq technology and assembled together with sequencing data sets obtained previously with Illumina GAIIx and 454 Roche FLX technologies ([@evu254-B72]). The assembly pipeline is summarized in [figure 1](#evu254-F1){ref-type="fig"}. F[ig]{.smallcaps}. 1.---Schematic representation summarizing the custom assembly pipeline used to reconstruct the genome of *Pseudomonas* sp. M1. Gray arrows indicate the pipeline stepwise flow. The differently colored lines associate each sequencing dataset with the phase(s) of the pipeline in which they were used.

MiSeq-sequencing yielded 5,431,166 paired-end reads of 250 bp. Raw reads were then corrected for sequencing errors, trimmed, quality-filtered, and approximately 70% of the original data set (3,847,594 reads) was de novo assembled using Edena ([@evu254-B34]) and Spades ([@evu254-B6]). Assembly-quality and number of SNVs and indels were assessed by mapping-filtered Illumina GAIIx paired-end data set (4,052,990 reads) against de novo draft genome assemblies using SNver ([@evu254-B82]). Merging Spades and Edena draft genomes with GAM-NGS assembly merger ([@evu254-B80]) and improvement with Pilon (<http://www.broadinstitute.org/software/pilon/>, last accessed October 1, 2014) yielded 52 contigs (longest contig size of 506,988 bp and N50 of 252,542 bp) with 9 SNVs and 12 indels. Scaffolding and gap closure were performed using the filtered data sets from Illumina MiSeq (3,847,594 paired-end reads) and GAIIx (4,052,990 paired-end reads and 4,698,623 mate-pair reads), which resulted in 20 contigs organized in two scaffolds with 9,303 unknown nucleotides (N's). After improvement using Pilon, mapping Sanger sequencing reads of regions evidencing abnormal coverage and/or gaps, mapping of the filtered 454 Roche FLX data set (253,532 single reads) and manual curation, the final improved genome draft of *Pseudomonas* sp. M1 was obtained in 4 contigs organized in 1 scaffold, with 21 unknown nucleotides (N's), 7 SNVs, 5 indels, and a total size of 6,982,449 bp. The resulting sequence was annotated with PGAP (Prokaryotic Genome Annotation Pipeline, NCBI) and deposited at DDBJ/EMBL/GenBank (accession number: ANIR00000000.2).

The improved draft of *Pseudomonas* sp. M1 genome includes 6,163 genes (6,023 coding sequences, 58 pseudogenes, 5 copies of each 5 S, 16 S, 23 S ribosomal RNA \[rRNA\] genes, 66 transfer RNA \[tRNA\] genes, and 1 noncoding RNA), with an average GC content of 67.2%. The major improvements of the current draft, when compared with the previous version, are reflected in the number of annotated genes and, importantly, in the reliability of base call along the genomic sequence, as deduced from the very low number of internal SNVs. Such improvement may have a significant impact on downstream functional analyses such as reference genome read mapping and protein identification.

The genome of *Pseudomonas* sp. M1 comprises a large set of oxygenases/hydroxylases-coding genes (99). Several enzymes belonging to this family have been shown to play an essential role on the initial oxidation steps of several natural and anthropogenic recalcitrant compounds. In silico prediction of M1 versatile metabolism by KAAS ([@evu254-B54]) suggested that the genomic repertoire of *Pseudomonas* sp. M1 may include the ability to engage in biotransformation/mineralization of recalcitrant compounds such as cresol and xylene, styrene, chloroalkanes and chloroalkenes, polycyclic aromatic compounds, naphtalene, benzoates, bisphenol, and atrazine ([supplementary table S1](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online). Additionally, M1 nutritional versatility is probably prompted by respiratory flexibility in using N-oxides as alternative electron donors, as suggested by the presence of denitrifying clusters in the genome. Such genomic features enhance the potential application of M1 for bioremediation at contaminated oxygen-depleted foci and as a biosensor for in situ monitorization of aliphatic and aromatic hydrocarbons.

Genome sequencing also allowed identification of several genetic features putatively involved in the metabolism of plant-derived volatiles other than β-myrcene ([supplementary table S1](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online), such as the acyclic terpene utilization cluster "atuRABCDEFGH" (*PM1_0225865--PM1_0225910*) that is described to be involved in citronellol and related terpenes catabolism in *Pseudomonas citronellolis* ([@evu254-B23]) and *Pseudomonas aeruginosa* strain PAO1 ([@evu254-B21]), and some oxygenases and dehydrogenases putatively involved in the metabolism of limonene- and alpha-pinene-related terpenes. Thus, *Pseudomonas* sp. M1 seems to harbor an attractive set of enzymes for biocatalysis of some high-value terpenoids.

RNA-Seq of M1 Cells Challenged with β-Myrcene
---------------------------------------------

Identification of genes involved in β-myrcene-stimulon was performed by comparative RNA-seq analysis of the transcriptomes of M1 cells grown in lactate (L) (control) and β-myrcene (M). Raw reads were trimmed and clipped for adaptor sequences and ambiguous nucleotides (N's) and mapped to a pool of rRNA and tRNA gene sequences from M1 draft genome to filter RNA-seq data sets exclusively for mRNA reads. To assess expression levels, mRNA reads were 99% uniquely mapped against M1 draft genome using EDGE-pro pipeline ([@evu254-B51]) and the Bioconductor package EdgeR ([@evu254-B63]). Among the large number of currently available tools for RNAS-seq analysis, this strategy was chosen based on: 1) EDGE-pro has been designed toward analysis of bacterial RNA-seq samples (in contrast with most of other read count applications which are directed toward the analysis of eukaryotic samples), and 2) EdgeR has been reported as one of the most used and reliable methods for data set normalization and statistical inference ([@evu254-B19]).

As shown in [table 1](#evu254-T1){ref-type="table"} (and more detailed in [supplementary table S2](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online), the comparative transcriptome analysis revealed a number of different functions being altered when cells are cultivated with β-myrcene. In M1 challenged cells, 627 genes were found to be upregulated at least 1.5-fold (all evidencing an FDR value lower than 0.05) and 1,246 genes were found to be downregulated at least 1.5-fold (all evidencing an FDR value of 0.05 or lower), suggesting an extensive molecular reprogramming of cell physiology. Globally, the cultivation of M1 cells in β-myrcene induced a change in the expression levels of over 30% of the genes in its genomic repertoire. Table 1Summary of Differential Expressed Genes in Myrcene-Grown *Pseudomonas* sp. M1 Cells, Organized in COG CategoriesFunctional CategoriesGenomeRNA-SeqNumber of GenesUp ≥ 1.5-FoldDown ≥ 1.5-Fold\[A\] RNA processing and modification110\[B\] Chromatin structure and dynamics502\[C\] Energy production and conversion35962108\[D\] Cell cycle control, cell division, chromosome partitioning3940\[E\] Amino acid transport and metabolism54333202\[F\] Nucleotide transport and metabolism100615\[G\] Carbohydrate transport and metabolism2111076\[H\] Coenzyme transport and metabolism1892635\[I\] Lipid transport and metabolism2435167\[J\] Translation, ribosomal structure and biogenesis1841215\[K\] Transcription4624463\[L\] Replication, recombination and repair14487\[M\] Cell wall/membrane/envelope biogenesis2252540\[N\] Cell motility178650\[O\] Posttranslational modification, protein turnover, chaperones1722119\[P\] Inorganic ion transport and metabolism2923590\[Q\] Secondary metabolites biosynthesis, transport and catabolism1601851\[R\] General function prediction only49040101\[S\] Function unknown4204645\[T\] Signal transduction mechanisms3012152\[U\] Intracellular trafficking, secretion, and vesicular transport1642241\[V\] Defense mechanisms621317\[W\] Extracellular structures101Total genes with COGs[^a^](#evu254-TF1){ref-type="table-fn"}4,463469974No COG assigned1,700158272Total6,1636271,246[^3]

Having in mind the assigned Clusters of Orthologous Groups (COG) categories, the most extensive alterations were registered in genes assigned to energy production and conversion; amino acid transport and metabolism; carbohydrate transport and metabolism; lipid transport and metabolism; transcription; inorganic ion transport and metabolism; and secondary metabolites biosynthesis, transport and catabolism. Importantly, β-myrcene metabolism seems to be dependent on expression of specific genetic loci, as deduced by a registered FC ratio greater than 10 for 52 of the upregulated genes.

*A Unique 28-kb Genomic Island Harbors the* myr+ *Genetic Trait*
----------------------------------------------------------------

β-myrcene supplementation induced the expression of 22 genes of a 28-kb genomic locus (*PM1_0216305--PM1_0216410*), from which 15 registered the highest expression levels from all β-myrcene-dependent transcriptome ([fig. 2](#evu254-F2){ref-type="fig"} and [supplementary table S2](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online). Interestingly, the *myrDABC* cluster (*PM1_0216385--PM1_0216400*), previously identified by [@evu254-B40] as essential for the initial steps of β-myrcene catabolism, is located in this 28-kb locus. Elements involved in conjugative transfer were found to be flanking this locus, which has an average GC content of 61.7% (genome average is 67.2%), indicating a hypothetical foreign nature for this genomic region. F[ig]{.smallcaps}. 2.---Global organization and functionality of the newly identified 28-kb GI harboring β-myrcene core-code. (*A*) Genomic locus organization and functional elements detected using in silico prediction. Genes transcribed in forward orientation are represented by rightward arrows, whereas genes located in antisense strand are represented by leftward arrows. Genes are organized in TU, promoter regions (P) are illustrated by blue arrows, and terminator sites indicated by red circles. (*B*) Detailed analysis of the RNA-seq data mapping with the GI. In the top panel, read mapping coverage per base was plotted for samples derived from M1 growth in lactate (gray line) and in β-myrcene (blue line). Prediction of the operon-like transcripts (horizontal orange bars) and transcription starting sites (TSS) (red arrows) was carried out with ReadXplorer ([@evu254-B36]). Expression levels for each gene are shown, in the bottom panel, as the ratio (FC) of normalized transcripts between cells grown in lactate (L) and cells grown in β-myrcene (M).

To characterize this GI, TU were defined in silico by predicting operon-like organization of genes from the GI, the respective promoter regions and terminator sequences. The 22 genes seem to be organized in eight TUs: TU1, TU2, TU3, TU4, TU5, TU6, TU7, and TU8 ([fig. 2](#evu254-F2){ref-type="fig"}*A*). Transcriptional organization was validated by mapping RNA-seq reads and de novo assembled contigs against the 28-kb GI, suggesting the polycistronic transcription of five of the eight TUs ([fig. 2](#evu254-F2){ref-type="fig"}*B*). Moreover, RNA-seq based analysis of this GI also allowed the prediction of seven TSS, six of which overlapping putative promoter regions deduced from in silico analysis ([fig. 2](#evu254-F2){ref-type="fig"}*A*). The detection of two TSS within TU3 may suggest that this TU may actually lead to different transcripts, depending on the growth conditions.

The first predicted TU comprises two acyl-CoA dehydrogenases (*PM1_0216305* and *PM1_0216310*) and an acetyl-CoA acetyltransferase*,* (*PM1_0216315*) putatively involved in β-myrcene catabolism, whose expression was increased approximately 365-, 370-, and 160-fold in β-myrcene-grown cells. A second cluster comprises two genes coding for a putative chemotaxis sensory transducer protein and a putative membrane protein (*PM1_0216320* and *PM1_0216325*). In the presence of β-myrcene, expression of this TU increased over 15-fold. From all 172 genes annotated on M1 genome as being involved in chemosensory systems and in signal transduction, 42 genes registered altered expression of at least 2-fold, from which only *PM1_0216320* and other signal protein-coding gene (*PM1_0226790*) were found to be upregulated in β-myrcene grown-cells, suggesting that TU2, together with *PM1_0226790*, might play essential role in recognition, transduction of the initial signal triggered by β-myrcene stimulus and possible uptake of β-myrcene (or derivative) molecules. OmpW family porins have been described to transport hydrophobic molecules across the membrane in some *Pseudomonas* spp. and *Escherichia coli* ([@evu254-B38]; [@evu254-B56]; [@evu254-B76]). Moreover, [@evu254-B56] showed that mutation in OmpW of *Pseudomonas fluorescens* resulted in loss of ability to catabolize polycyclic aromatic molecules, such as naphtalene. TU3 codes for six products that might also be involved in β-myrcene oxidation: Another acetyl-CoA acetyltransferase (*PM1_0216330*) with increased expression of approximately 320-fold, a monooxygenase (*PM1_0216335*) increased by 431-fold, an enoyl-CoA hydratase (*PM1_0216340*) increased by 412-fold, an acyl-CoA dehydrogenase family protein (*PM1_0216345*) upregulated by 263-fold, an epoxide hydrolase (*PM1_0216350*) overexpressed by 102-fold, and a hypothetical protein (*PM1_0216355*) with increased expression by 51-fold. Moreover, TU4 is composed by a single gene, coding for a putative LuxR-family transcriptional regulator (*PM1_0216360*). This putative regulator registered increased expression of 190-fold in β-myrcene-grown cells. The fifth transcriptional unit is composed by two genes which code a rubredoxin (*PM1_0216365*), whose expression was increased by 260-fold, and a fatty acid desaturase (*PM1_0216370*), whose expression was increased by 202-fold. The sixth TU, comprises five genes which code for another LuxR-family transcriptional regulator (*PM1_0216375*), an oxidoreductase (*PM1_0216380*), the acyl-CoA synthase *myrC* (*PM1_0216385*), the alcohol dehydrogenase *myrB* (*PM1_0216390*), and the aldehyde dehydrogenase *myrA* (*PM1_0216395*). Expression levels were upregulated by 5-fold, 71-fold, 129-fold, 86-fold, and 44-fold, respectively. Gene coding for the enoyl-CoA hydratase *myrD* (*PM1_0216400*) comprises a seventh transcriptional unit and registered an increased expression of 220-fold. The eighth transcriptional unit predicted is composed by two genes which code for a putative LysR-family transcriptional regulator (*PM1_0216405*) and a hypothetical protein (*PM1_0216410*), whose expression was increased by 95- and 70-fold, respectively.

The high expression of such set of catabolic enzymes strongly suggests that this GI is the main genetic trait that provides *myr +*phenotype to *Pseudomonas* sp. M1. Particularly, the monooxygenase-coding gene *PM1_0216335* showed the highest expression registered in all transcriptome and might be involved in initial oxidation of β-myrcene into a more soluble substrate. The understanding of its substrate range would be of utmost interest for terpene-based industry and future fine-tuning of its activity for biotechnological applications. Moreover, the differential expressed genes encoding for putative regulators (*PM1_0216360, PM1_0216405*, and *PM1_0216375*) suggest that they may be key players in the modulation of expression of the catabolic enzymes encoded within the GI.

Expression of predicted TUs was validated by RT-qPCR according to Pfaffl ([@evu254-B58]) and using *recA* (*PM1_02298*) as reference gene for normalization ([@evu254-B52]; [@evu254-B29]) ([fig. 3](#evu254-F3){ref-type="fig"}). Similar to expression levels measured by RNA-seq, key genes *PM1_0216305, PM1_0216325*, *PM1_0216335*, *PM1_0216350*, *myrA,* and *myrD* were highly expressed in M1 cells challenged with β-myrcene, as showed by the FC ratio in [figure 3](#evu254-F3){ref-type="fig"}. In general, the FC ratios obtained by RT-qPCR were higher when compared with RNA-seq-based estimated FC ratios. Although both methods have been proved to provide significant information regarding the evaluation of mRNA levels among different conditions, they rely on different methodological approaches which can differently affect the estimation of FC ratios (e.g., library preparation, data normalization, amount of starting template, and template amplification procedures). Nonetheless, both approaches clearly show that the expression of the selected genes is strongly stimulated by β-myrcene. F[ig]{.smallcaps}. 3.---Expression levels of key genes from the 28-kb GI detected by RNA-seq and validated by RT-qPCR, shown as the ratio of normalized transcripts between cells grown in lactate and cells grown in β-myrcene (FC).

The 28-kb GI Evolved from the Assembly of Functional Blocks
-----------------------------------------------------------

As the GC content percentage within this GI has an average of 61.7% which is significantly different of the whole-genome average (67.2%), it was hypothesized that this GI may have been acquired as a whole DNA sequence segment through horizontal gene transfer. However, a thorough inspection of the public databases for DNA sequence deposition indicated that this GI was not present in any records besides the ones related with *Pseudomonas* sp. M1. In fact, the closest homologs of each predicted gene of the GI indicated a wide range of percentage of identity variation ranging from 56% to 87% (for the top hits), *PM1_261405* evidencing a closest similarity with the matched homologs (87%) and *PM1_216375* evidencing a lowest similarity (56%). To obtain clues on the origin of this GI, a phylogenetic analysis of all the genes of the 28-kb GI identified in this work was carried out. In this context, and due to the registered fairly low identity percentage between the genes of the GI and the respective homologs, only the 20 closest related sequences for all the genes (except *PM1_216355* for which only ten sequences resembling this gene were found) were selected for this analysis. Additionally, the analysis of the *guaA* gene of M1 strain and the respective 20 closest sequences was also performed. This gene was selected as a representative gene of M1 core genome as it evidences a GC content % similar to the average of the whole genome of *Pseudomonas* sp. M1 and it often used in multilocus sequence typing of *Pseudomonas* strains. In [figure 4](#evu254-F4){ref-type="fig"}, and in more detail in [supplementary table S3](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online, the results of the comparison of each gene of GI with their homologs are presented. The comparison of each gene of the GI with the respective 20 homologs indicates 1) a significantly lower average percentage of identity (57--77%) with respect to the 93% of identity for the reference gene *guaA* and 2) a significant difference in the GC% content of M1 genes with respect with the closest homologs (ranging 1--12% difference). These differences are supported, at least partially, by the diversified taxonomic origin of the closest homolog sequences, as depicted in [figure 4](#evu254-F4){ref-type="fig"}. As expected, the results for the representative *guaA* gene are in agreement with previously published phylogenetic analysis of *Pseudomonas* sp. M1 based on 16 S rRNA sequencing ([@evu254-B67]) and groups M1 within the Gammaproteobacteria class (and *Pseudomonas* genus). In contrast, most of the genes of the GI, with the exception of *PM1_216320*, share significant identity with homolog sequences from organisms belonging to a taxonomic class not related with *Pseudomonas* sp. M1. In the case of *PM1_216305*, *PM1_216310*, *PM1_216315*, *PM1_216320*, *PM1_216325*, *PM1_216330*, *PM1_216335*, *PM1_216340*, *PM1_216355*, *PM1_216360*, *PM1_216365*, *PM1_216375*, *PM1_216405*, and *PM1_216410* genes, 75% (or more) of their homolog sequences belong to organisms of a specific taxonomic class. Distinctively, in the case of *PM1_216345*, *PM1_216350*, *PM1_216370*, *PM1_216380*, *PM1_216385*, *PM1_216390*, *PM1_216395,* and *PM1_216400* genes, the taxonomic distribution of their homolog sequences is spread among different classes. This asymmetric taxonomic distribution of the homolog sequences of the GI may suggest that the genetic reservoirs of some of the genes within the GI are confined to specific taxonomic groups whereas homologs of the other genes are widely distributed across different taxonomic classes. Moreover, this asymmetry also suggests that, probably, the GI resulted from acquisition of DNA segments from different environmental reservoirs which were assembled together evolving toward gene specialization for β-myrcene utilization. This hypothesis is in agreement with the multiple parallel evolutionary origins of GIs outlined before ([@evu254-B41]). F[ig]{.smallcaps}. 4.---Graphical representation of the analysis of the results of multialignments of each gene of the GI implicated in β-myrcene catabolism in *Pseudomonas* sp. M1 and PM1_216145 (*guaA*), and the respective most homolog genes, retrieved using BLASTn. The dotted line indicates the average percentage of identity with respect to M1 genes. The average of GC percentage of each M1 gene or of the homologs sequences is represented by blue or red lines, respectively. The taxonomic origin of each homolog sequence was assessed and grouped according the respective taxonomic class. The bars in graph represent relative distribution (in percentage of the total number homolog sequences) of Alphaproteobacteria (green), Actinobacteria (light blue), Betaproteobacteria (red), Gammaproteobacteria (dark blue), and other classes (light orange).

Each gene of the GI, the *guaA,* and the most the selected most homolog sequences (analyzed above) were aligned using MAFFT. The resulting multialignments were used to construct a phylogenetic tree for each gene (GI and *guaA*) using an ML method. The respective −Log likelihood (−ln *L*) scores and divergence distance with respect to each *Pseudomonas* sp. M1 gene (estimated as the average of the patristic distances between tree tips and the respective M1 gene) were extracted. The −ln *L* score reflects the overall phylogenetic tree topology and branch lengths as a consequence of the relatedness of the sequences under analysis, being higher values an indication of higher phylogenetic distances. In fact, in this study it was verified that, in general, the −ln *L* increased when the number of taxonomic genus represented in the multialignment sequences increased (data not shown). Presumably, M1 genes evidencing closer −ln *L* scores and divergence distances, based on the comparison with the selected homolog sequences, may share more similar evolutionary histories. Therefore, a multivariable matrix including the −ln *L* scores and divergence distances for all the genes of the GI and *guaA* was constructed. In [figure 5](#evu254-F5){ref-type="fig"}, the agglomerative clustering analysis of the resulting matrix is shown. As expected from the results presented in [figure 4](#evu254-F4){ref-type="fig"}, the reference gene *guaA* clustered distantly from the genes within the GI as the homolog sequences share a higher degree of identity. Remarkably, this analysis suggests that most of the genes within TU3 and TU6 ([fig. 2](#evu254-F2){ref-type="fig"}) may share a close evolutionary history as suggested by their inclusion in the same cluster partition, C9 and C6 of [figure 5](#evu254-F5){ref-type="fig"}, respectively. Additionally, *PM1_216320/PM1_216325* and *PM1_216405/PM1_216410* pairs of genes cluster close to each other and correspond to the predicted TU2 and TU8, respectively. These observations are in agreement with the correlation between coexpression and coevolution of genes encoding enzymes with coupled metabolic fluxes previously reported (Seshasayee et al. 2009). Moreover, *PM1_216320, PM1_216360,* and *PM1_216375* genes, that putatively code for three of the four transcriptional regulators present in the GI, cluster closer (clusters C1, C2, and C3 of [fig. 5](#evu254-F5){ref-type="fig"}, respectively) and apart from the catabolic genes. Significantly, the highest values for the −ln *L* and/or average patristic distance parameters were registered for these genes. Our results may indicate that the control of the catabolic genes of the GI requires an intricate tuning of the regulatory sequences towards an effective responsiveness to β-myrcene (and possibly other terpenes), probably due to the unusual nature of the substrates. Thus, it suggests that the regulatory and the catabolic components of the GI may have evolved separately. This hypothesis is in agreement with a modular assembly concept of emerging catabolic pathways that includes evolution from constitutive to regulated expression of the catabolic enzymes, as previously described ([@evu254-B14]). F[ig]{.smallcaps}. 5.---Multivariate clustering analysis of the inferred phylogenetic tree-based parameters −log likelihood score (−ln *L*) and average patristic distance (Divergence) for each gene of the GI. For a comparison purpose, the parameters for the representative gene *guaA* were also included in the analysis. Each M1 gene, represented in the figure, was multialigned together with the most homolog sequences and clustering analysis of the inferred parameters was carried by estimating the squared Euclidean distance between the different sequences in each multialignment, followed by agglomerative clustering with UPGMA method. The number of clusters was estimated by partitioning around mediods clustering algorithm. Each cluster is outlined by a different color.

β-Myrcene-Induced Genetic Repertoire Might Code for Enzymatic Redundancy toward a Broader Range of Terpene-Substrates
---------------------------------------------------------------------------------------------------------------------

Transcriptome data not only support the pathway described by [@evu254-B40] but also strongly suggest the involvement of several other gene products in β-myrcene oxidation and raise new hypotheses about β-myrcene metabolic pathway that will be the focus of following studies. At this moment, it is unclear whether all 22 genes from the 28-kb GI are actually necessary for the metabolism of the aliphatic C8-chain of β-myrcene or whether the high expression levels of the island are induced by a chemosensory signal common to a broader metabolic potential toward several terpene-backbone compounds. The expression of the putative epoxide hydrolase (*PM1_0216350*) when *Pseudomonas* sp. M1 is grown in β-myrcene medium may indicate a possible epoxidation of β-myrcene for subsequent hydrolysis by PM1_0216350 into a β-myrcene alcohol. However, terpene-epoxide microbial biotransformation has been reported mainly as result of monooxygenase activity toward cyclic molecules (e.g., alpha-pinene and limonene) ([@evu254-B79]). Therefore, it is also possible that the putative epoxide hydrolase encoded by the GI may be actually involved in epoxidation of substrates other than β-myrcene.

β-myrcene biotransformation has been reported in other gram-negative bacteria. *Rhodococcus erythropolis* MLT1 was described to biotransform aerobically β-myrcene into geraniol ([@evu254-B75]). In *Castellaniella defragans* 65Phen β-myrcene is hydratated to linalool, following isomerization to geraniol and oxidation to geranic acid, under anaerobic conditions ([@evu254-B12]; [@evu254-B48]; [@evu254-B47]; [@evu254-B49]), which indicates a cross-talk between catabolic pathways of different terpenes. In M1 cells challenged with β-myrcene, the majority of genes from the acyclic terpene utilization cluster *atuRABCDEFGH* and the leucine/isovalerate degradation cluster *liuRABCDE* (*PM1_0223355*--*PM1_0223380*) were found to be upregulated (average cluster expression induced by 3.8-and 4.1-fold, respectively). The overexpression of these two clusters suggests that, similarly to *R. erythropolis* and *C. defragans* ([@evu254-B21]; [@evu254-B23]), a citronellol/geraniol-like compound might be an intermediate/byproduct of β-myrcene metabolism and channel β-myrcene-derivatives through propionyl-CoA into central metabolism ([fig. 6](#evu254-F6){ref-type="fig"}). Another hypothesis is that the sensory system present in the GI may be involved in the recognition of a spectrum of terpene-backbone molecules, including β-myrcene, and trigger expression of several genetic loci with complementary/redundant degradative functions ([fig. 6](#evu254-F6){ref-type="fig"}). F[ig]{.smallcaps}. 6.---Schematic prediction of β-myrcene metabolic pathway and redistribution of central metabolic fluxes in M1 cells challenged with β-myrcene. Upregulated genes are shown by upward black (\>2-fold) or gray (\<2-fold) arrows and downregulated genes are shown by downward black (\>2-fold) or gray (\<2-fold) arrows. Genes with similar expression levels in *Pseudomonas* sp. M1 cells grown in lactate medium and β-myrcene medium are indicated by white squares.

Several other genes located outside of the 28-kb locus and coding for beta-oxidation-like enzymes had their expression altered in M1 cells grown in β-myrcene medium. As shown in [supplementary table S2](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online, increased expression by at least 2-fold was detected for 8 acyl-CoA-like dehydrogenases (*PM1_0206520, PM1_0208500, PM1_0211935, PM1_0220850, PM1_0221270, PM1_0224730, PM1_0224740,* and *PM1_0224920*), a putative acetyl-CoA acetyltransferase (*PM1_0224735*), a putative acetyl-CoA synthetase (*PM1_0208785*), 2 putative acyl-CoA thiolases (*PM1_0221265* and *PM1_0212440*), 2 putative enoyl-CoA hydratases (*PM1_0221070* and *PM1_0221075*), a putative enoyl-CoA reductase (*PM1_0225900*), and 2 putative aldehyde dehydrogenases (*PM1_0226785* and *PM1_0230170*). In β-myrcene-grown *C. defragans* and *R. erythropolis*, overexpression of additional acyl-CoA dehydrogenases and oxidoreductases was detected ([@evu254-B75]; [@evu254-B49]). Furthermore, kinetic assays showed that geraniol dehydrogenase GeoA of *C. defragans* was able to catalyze either cyclic (e.g., perillyl-alcohol) and acyclic monoterpenes (e.g., geraniol) ([@evu254-B47]), and its deletion did not impair completely growth on β-myrcene medium, suggesting the existence of secondary enzyme system acting on acyclic terpene-backbone molecules ([@evu254-B47]). Similarly, mutagenesis trials in *Pseudomonas* sp. M1 (e.g., *ΔmyrB* mutant) did not abolish completely the ability of M1 to grow in β-myrcene medium (data not shown), supporting the redundancy of beta-oxidation-like enzymes overexpressed during growth in β-myrcene medium and suggesting a cross-talk between metabolic fluxes of terpene catabolism and a complementary substrate range for these set of enzymes ([fig. 6](#evu254-F6){ref-type="fig"}).

Central Metabolic Flux Is Changed in β-Myrcene-Grown Cells
----------------------------------------------------------

Concomitant with the upregulation of the 28-kb, expression of almost all genes coding enzymes of the Entner--Doudoroff pathway was downregulated ([supplementary table S2](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online). On the other hand, genes coding for enzymes of the tricarboxylic acid cycle (TCA) show similar expression levels in M1 cells grown in lactate or β-myrcene ([fig. 6](#evu254-F6){ref-type="fig"} and [supplementary table S2](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online), with exception of the citrate conversion to isocitrate as the *acnA* gene (*PM1_0223060*) coding for aconitate hydratase registered 2.4-fold increase in expression. Subsequently, isocitrate might be channeled to the glyoxylate shunt, mediated by a 2.6-fold increase in expression for *PM1_0221765* (coding for a isocitrate lyase AceA-like). The glyoxylate bypass might increase levels of succinate for energy production, surpassing the metabolic branching of TCA intermediaries (e.g., 2-oxoglutarate) for amino acid biosynthesis. Also, glyoxylate shunt might increase levels of malate, which is further converted to oxaloacetate, a key branching step in the central metabolism. *Pseudomonas* sp. M1 genome harbors *mqoB* gene (*PM1_0214580*) coding for a malate:quinone oxidoreductase and a gene coding for a putative malate dehydrogenase (*PM1_0229630*), being the former upregulated 2.0-fold in M1 cells grown in β-myrcene. In *P. citronellolis*, *mqoB* expression has been described to be essential for the metabolism of the acyclic terpenes citronellol and citronellic acid ([@evu254-B22]). A possible bypass of isocitrate toward the glyoxylate shunt in β-myrcene-grown cells supports previous proteomic data ([@evu254-B68]), which suggested that TCA flux might be directed toward oxaloacetate synthesis as this intermediate might be used in condensation reaction with a β-myrcene derivative, such as propionyl-CoA ([@evu254-B74]). This would channel the final intermediates of β-myrcene catabolism to the central metabolism through 2-methycitrate cycle (e.g., in form of succinate and pyruvate) ([fig. 6](#evu254-F6){ref-type="fig"}). In fact, the expression of all five genes coding for 2-methycitrate cycle enzymes (*PM1_0215080--PM1_0215100*) was also increased in β-myrcene-grown cells (average upregulation of 10.6-fold).

β-Myrcene Catabolism Induces an Anaerobiose-Like Physiology
-----------------------------------------------------------

RNA-seq data also suggested that M1 cells may use alternative pathways for energy production during β-myrcene catabolism, resembling an oxygen-depletion-like physiology ([supplementary table S2](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online). Expression levels of four clustered genes coding for Dnr-like regulators (*PM1_0227665, PM1_0227670, PM1_0227660,* and *PM1_0227675*) were upregulated 7.7-, 2.5-, 2.1-, and 1.9-fold, respectively. Furthermore, Dnr-like regulators were probably responsible for inducing the transcription of several denitrification operons: Average 2.3-fold induction of the *narK1K2GHJI* (*PM1_0215670--PM1_0215695*), an average 38.5-fold induction of the *norCBD* cluster (*PM1_0227680--PM1_0227690*), an average induction of 6.2-fold for the *nirQOP* cluster (*PM1_0227695--PM1_0227705*), 7.1-fold increase in expression of *nirM* (*PM1_0227710*), and an average induction of 3.3-fold for *nosRZDFYL* cluster (*PM1_0227715--PM1_0227740*), coding the set of reductases and cytochromes necessary for N-oxides respiration ([@evu254-B4]). Moreover, expression of genes from the Anr-dependent operon *arcDABC* (*PM1_0202140--PM1_0202155*), coding for enzymes of the arginine deiminase pathway, was also induced at least 4-fold in β-myrcene-grown M1 cells. This operon is responsible for arginine fermentation with the concomitant formation of ATP ([@evu254-B27]). Arginine fermentation might have been induced at earlier stages as a source of energy production that would be gradually replaced by nitrate respiration as Anr--Dnr transcriptional cascade activated their target genes and as β-myrcene metabolism demanded higher yields of energy production. Catabolism of some terpenes in nitrate-respiring conditions has already been reported for *C. defragrans* (e.g., alpha-pinene, β-myrcene, and limonene) ([@evu254-B25]; [@evu254-B35]; [@evu254-B12]) and *Thauera* spp. (e.g., linalool, menthol, and eucalyptol) ([@evu254-B24]), mainly based on metabolite identification. In this study, a molecular snapshot shed some light in what might be a common mechanism for terpene-backbone metabolism under denitrifying conditions.

M1 Cells Reprogram Their Physiology during β-Myrcene Catabolism
---------------------------------------------------------------

As described above, β-myrcene stimulus activates specific pathways on *Pseudomonas* sp. M1, which in turn appears to switch off genes that might become a metabolic burden for growing cells. M1 cells grown in β-myrcene show several secondary metabolic pathways and respective transport systems downregulated when compared with cells grown in lactate ([supplementary table S2](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online), mainly those related to amino acids biosynthesis (e.g., glutamate, glutamine, serine, histidine, lysine, arginine, proline), utilization of alternative carbon sources (e.g., aromatic compounds, sugars), degradation of urea, amines, and polyamines.

The altered expression of genes involved in lipopolysaccharide (LPS) and exopolysaccharide synthesis suggested a reorganization of cell envelope in β-myrcene grown cells. Expression of a cluster harboring the *wzm-wzt-wbpXYZ* genes (*PM1_0228115--PM1_0228170*) ([@evu254-B64]; [@evu254-B28]) involved in LPS biogenesis was downregulated over 1.5-fold, as well as a gene coding for the LptA protein (*PM1_0228970*), an element of the LPS transport system across the outer membrane ([@evu254-B77]), whose expression was downregulated by 1.4-fold. The expression of the main cluster associated with alginate synthesis, *algD-8-44-KEGXLIJFA* (*PM1_0215545--PM1_0215600*) ([@evu254-B15]) was downregulated about 2-fold, suggesting a reduction in the synthesis of this polymer. Changes in outer cell envelope were, apparently, followed by changes in membrane lipid composition, membrane stability and permeability as many genes coding for the synthesis of lipid polar head group, transporters/facilitators, porins and membrane proteins had their expression changed when comparing with lactate grown-cells. In β-myrcene-grown cells, expression of two genes coding for OprD homologous (*PM1_0203585* and *PM1_0218170*) was increased by 70.8-fold and decreased by 10.0-fold, respectively. Moreover, expression of a set of efflux pumps was increased more than 2-fold, including six genes coding for Resistance-Nodulation-Cell Division family efflux systems (*PM1_0222160, PM1_0210745, PM1_0222165, PM1_0217125, PM1_0217120,* and *PM1_0202020*), a gene homologous to fusaric acid-inducible resistance pumps (*PM1_0210750*) and two genes coding for a multidrug resistance transporter system (*PM1_0202030* and *PM1_0202025*) with homology to *emrAB* cluster of *E. coli*. On the other hand, OmpQ-family protein-coding genes *PM1_0230840* and *PM1_0213495* were downregulated to 1.8- and 2.0-fold, respectively; *opdE* (*PM1_0220290*) was downregulated to 1.8-fold; and *opdH* (*PM1_0216660*) expression was reduced to 2.0-fold.

Abundance of porins changes in response to nutritional requirements of the cell and control membrane permeability toward stressors ([@evu254-B84]; [@evu254-B73]; [@evu254-B62]). Additionally, efflux systems are usually described to be upregulated in stressful conditions, being responsible for cell homeostasis by exporting reactive and toxic compounds ([@evu254-B26]; [@evu254-B60]; [@evu254-B39]). Furthermore, such cell envelope dynamics might account for the adaptive mechanisms toward β-myrcene hydrophobicity, which might affect overall membrane integrity and induce cellular stress ([@evu254-B71]; [@evu254-B78]; [@evu254-B16]). Increased membrane hydrophobicity was previously described in *Pseudomonas* strains challenged by other hydrophobic compounds, such as n-alkanols and phenanthrene, as a mechanism to enhance biotransformation ([@evu254-B32]; [@evu254-B7]). Overall, M1 cells grown in β-myrcene, apparently, leads to a reorganization of membrane structure and functions, most probably to cope with the highly hydrophobic character of β-myrcene and eventual generation of β-myrcene-derived toxic intermediates, as previously evidenced ([@evu254-B68]). Additionally, changes in M1 cell envelope structure and functions may account for a strategy to increase substrate availability/contact surface at the oil--water interface by changing cell surface polarity and control levels of stress by changing membrane permeability and stiffness.

Detailed characterization of expression kinetics and regulation are required to fully understand the β-myrcene-stimulon. The sensory system of β-myrcene shows high sensitivity in detecting β-myrcene highly hydrophobic molecules, being an interesting candidate for biosensory applications. Currently, three enzymes, the monooxygenase PM1_0216335, the epoxide hydrolase PM1_0216350 and the 2-methylcitrate synthase PrpC (PM1_0215095), are in protein crystallization trials (data not shown) for future fine-tuning and source of novel biomolecules. Nevertheless, key players involved in the initial steps of β-myrcene metabolism by *Pseudomonas* sp. M1 and the range of plant-derived volatiles that may be biotransformed by the 28-kb GI, as a functional block and by M1 as a cell-factory, still remain to be elucidated and will be focus of following work.

Supplementary Material
======================

[Supplementary tables S1--S3](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) are available at *Genome Biology and Evolution* online (<http://www.gbe.oxfordjournals.org/>).

###### Supplementary Data

This work was supported by FEDER through POFC---COMPETE and by national funds from Foundation for Science and Technology (Portugal) through the projects PEst-C/BIA/UI4050/2011, PTDC/EBB-BIO/104980/2008 and PTDC/BIA-MIC/113733/09, and through a PhD grant (grant number SFRH/BD/76894/2011) to P.S.-C.

[^1]: **Data deposition:** The Whole Genome Shotgun project of *Pseudomonas* sp. M1 has been deposited at DNA Data Bank of Japan/EMBL/GenBank under the accession ANIR00000000. The version described in this paper is ANIR00000000.2.

[^2]: **Associate editor:** José Pereira-Leal

[^3]: ^a^In total, 482 genes of the whole-genome annotation have multi-COG assignment. Further details are provided in [supplementary table S2](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1), [Supplementary Material](http://gbe.oxfordjournals.org/lookup/suppl/doi:10.1093/gbe/evu254/-/DC1) online.
