Transcriptomic Analyses Unveil Hydrocarbon Degradation Mechanisms in a Novel Polar Rhodococcus sp. Strain R1B_2T From a High Arctic Intertidal Zone Exposed to Ultra-Low Sulphur Fuel Oil

Abstract

As Arctic shipping increases due to climate change, characterised by rising temperatures and decreasing sea-ice coverage, the risk of oil spills through the Northwest Passage in this fragile ecosystem grows, necessitating effective bioremediation strategies. Research on bioremediation using Arctic coastal sediment bacteria has gained attention, particularly Rhodococcus species that play key roles in hydrocarbon degradation under extreme conditions. This study investigates the hydrocarbon degradation capabilities of a novel cryophilic Arctic Rhodococcus sp. strain R1B_2T isolated from Canadian high Arctic beach sediment in Resolute Bay, exposed to ultra-low sulphur fuel oil for 3 months at 5°C. Comparative transcriptomics analyses revealed dynamic responses and metabolic plasticity, with upregulation of genes for aliphatic, aromatic, and polycyclic aromatic hydrocarbons, biosurfactant production (e.g., rhamnolipid), cold adaptation, and stress responses. The strain possesses several key alkane degradation genes (alkB, almA, CYP153, ladA), with co-expression network analysis highlighting synergistic mechanisms between alkB and CYP153 that target different chain-length alkanes (alkB: ~C5-C20; CYP153: ~C5-C12 and > C30), demonstrating complementary degradation strategies. The findings reveal adaptive mechanisms and degradation kinetics of native Arctic bacteria, highlighting the potential of Arctic cryophilic and halotolerant Rhodococcus species for oil spill remediation in polar marine environments.

Keywords: Rhodococcus; Arctic bacteria; biodegradation; comparative transcriptomic; differential gene expression; hydrocarbon degradation.

FIGURE 1

Novel strain R1B_2T phylogenetic position compared to 55 other species of Rhodococcus. (A) A Maximum Likelihood tree generated using RAxML with 16S rRNA gene sequences of 55 strains and species of Rhodococcus was aligned along with one reference sequence of a short length 16S rRNA gene (sequenced using the Sanger method). The phylogenetic tree was constructed from an alignment of 56 sequences, spanning 909 bp, with bootstrap support calculated over 1000 repetitions. Only bootstrap values greater than 50 (out of 100) are displayed. (B) Heatmap of average nucleotide identity (ANI) between 8 referenced genomes of Rhodococcus species and novel strain R1B_2T. ANI and genome information can be found in Table S3.

FIGURE 2

Petroleum hydrocarbon analyses of Rhodococcus exposed to ultra‐low sulphur fuel oil (ULSFO) over 3 months. (A) and (B) Alkane fractions (F2 to F4) degradation over time. (C) and (D) Aromatic compounds degradation over time. Error bars are for the triplicate cultures incubated for 1 and 3 months with biological replicates (n = 6). Missing error bars mean that triplicate values were not available. A complete list of values is in Table S1. Details of the One‐way ANOVA test for natural attenuation by Rhodococcus sp. strain R1B_2T over time are in Table S5.

FIGURE 3

Clustering of transcriptome samples. (A) Hierarchical clustering of all samples based on normalised reads and unweighted pair group method with arithmetic mean (UPGMA) based on the Bray‐Curtis distance matrices of significant DEGs for all sampling times. (B) Canonical correspondence analysis (CCA) with the Bray–Curtis dissimilarity measure of reads per gene of all transcriptomes. (C) Top 25 differentially expressed genes of the Rhodococcus sp. strain R1B_2T. Heatmap of clustered results based on normalised and transformed read counts. The functional annotation of each gene is on the right of the heatmap, indicated letters a to c, refer to Table S6. The variance stabilising transformation (vst) is calculated with the vst() function within the DeSeq2 package on R and transforms the count data to normalise the read count. The size factor corresponds to the sequencing depths of each transcriptome.

FIGURE 4

Top expression of differentially expressed genes (DEGs) annotated with gene ontology (GO) terms between comparisons of time and with or without ULSFO. (A) GO term of the biological process category. (B) GO term of the Cellular Component category. (C) GO terms of the Molecular Function category. A complete list of genes and annotations can be found in Table S7.

FIGURE 5

Weighted gene co‐expression network analysis (WGCNA) of differentially expressed genes (DEGs). (A) The dendrogram shows the DEGs that clustered into 8 modules. Genes were clustered based on a dissimilarity measure (1‐TOM). (B) Module with weighted correlations and corresponding p values at time of sampling, and salinity conditions. The colour scale shows module‐trait correlation based on Pearson's rank correlation. Asterisks within the correlation indicate the significance level p value from Pearson's Rho (p value < 0.05 *, < 0.01 ** and < 0.001 ***). (C) Scatter plots showing correlation between module membership and gene significance for four key modules (brown, blue, turquoise and grey). Each point represents a gene; trend lines show linear relationship; Cor = correlation coefficient, p = significance value. (D) Brown, turquoise and blue modules co‐expression network based on WGCNA analysis representing the significant DEGs interactions. The node size is calculated on the edge count. The higher value or bigger size node refers to a stronger connection or co‐expression of DEGs. Black dotted circles represent the hub gene in each module. Group of nodes from brown module are referred as A1 to A5, from blue module as B1 to B5, from turquoise module as C1 to C5 and from grey module as D1. Complete list of gene in each node and module can be found in Table S10.

FIGURE 6

Circular heatmap of differentially expressed genes (DEGs) and their expression profile based on the log₂ fold change between time of sampling (T0, T1—1 month, T3—3 months) and with or without ULSFO. Colours in the outer circle indicate alkane or aromatic compound degradation category for each DEG. The heatmap uses colour intensity to show the level of gene expression, with darker colours indicating higher expression levels. A complete list of genes and annotations can be found in Table S11.

FIGURE 7

Key genes expressed in Rhodococcus for alkane compound degradation. A complete list of genes and annotations can be found in Tables S11 and S13.

FIGURE 8

Key genes expressed in Rhodococcus for aromatic compound degradation. A complete list of genes and annotations can be found in Table S11.