Uncovering the temporal dynamics and regulatory networks of thermal stress response in a hyperthermophile using transcriptomics and proteomics

Abstract

Extreme environments provide unique challenges for life, and the study of extremophiles can shed light on the mechanisms of adaptation to such conditions. Pyrococcus furiosus, a hyperthermophilic archaeon, is a model organism for studying thermal stress response mechanisms. In this study, we used an integrated analysis of RNA-sequencing and mass spectrometry data to investigate the transcriptomic and proteomic responses of P. furiosus to heat and cold shock stress and recovery. Our results reveal the rapid and dynamic changes in gene and protein expression patterns associated with these stress responses, as well as the coordinated regulation of different gene sets in response to different stressors. These findings provide valuable insights into the molecular adaptations that facilitate life in extreme environments and advance our understanding of stress response mechanisms in hyperthermophilic archaea.

Keywords: archaea; cold shock; heat shock; proteomics; transcriptomics.

Fig 1

Validation of the experimental setup by temporal analysis of the heat and cold shock response in Pyrococcus furiosus via transcriptomics and proteomics. (A) Experimental design illustrating the seven stress-related conditions, including cold shock (CS) at 4°C and heat shock (HS) at 105°C, along with recovery samples, analyzed using RNA-seq and proteomics, compared to the control condition (Ctrl) at the optimal growth temperature of 95°C. (B) Venn diagram displaying the overlap between identified protein-coding transcripts (gray) and total identified proteins (yellow) under stress-related and control conditions. Genes only identified using RNA-seq or proteomics are shown at the left and right side, respectively. (C) Distribution of transcripts per million (TPM) normalized transcript counts, color-coded based on the legend in panel A. (D) Distribution of intensity-based absolute quantification (iBAQ) values for proteomics samples. (E) Principal component analysis of RNA-sequencing samples based on total counts; replicates are indicated using different symbols, with outliers removed. (F) Circular genome-wide plot of the P. furiosus genome (position 0 at top) with protein-coding genes on the two outer rings; Z-score normalized TPM counts color-coded by dark brown (negative) to dark green (positive) for each condition in the following order (outer to inner): Ctrl, CS 1, CS 2, CS 3, CS R, HS 1, HS 2, HS R. (G) PCA and (H) genome-wide circular plot for mass spectrometry samples analogous to panels E and F.

Fig 2

Immediate heat shock leads to extensive reprogramming on the transcriptome level with moderate overlap to cold shock response. (A) Volcano plot displaying log₂-fold changes (x-axis) and significance (−log₁₀ adjusted P-value, y-axis) comparing transcriptome changes at HS 1 (5 min) to the control condition. Protein-coding transcripts are categorized by significance and fold changes, color-coded as strongly upregulated (padj < 0.05 and log₂FC ≥ 1, dark green, bold font), upregulated (padj < 0.05 and log₂FC < 1 and log₂FC > 0, light green, normal font), non-regulated genes (NS, padj ≥ 0.05, white), strongly downregulated (padj < 0.05 and log₂FC ≤ −1, dark brown, bold font), downregulated (padj < 0.05 and log₂FC > −1 and log₂FC < 0, light brown, normal font). HS genes are highlighted. Significance level of 0.05 is shown as a dashed line. (B) Schematic representation of the HS response in P. furiosus, with genes color-coded according to transcriptome fold changes. (C) Gene set enrichment analysis of archaeal clusters of orthologous genes (arCOGs), with significance levels indicated by a color bar ranging from white to dark purple. Genes with a P-value <0.05 are considered significantly overrepresented and highlighted by a circle, with circle size reflecting the number of genes in the category. (D) Volcano plot showing log₂-fold changes for CS 1 condition compared to the control condition, with the y-axis displaying significance. A cutoff value of adjusted P-value of 0.05 is indicated by a dashed line, and relevant genes are highlighted. (E) Functional enrichment analysis using arCOG descriptions; significance levels and the number of regulated genes are shown in the upper legend. (F) Overlap between strongly upregulated (dark green) and downregulated (dark brown) genes in HS 1 and CS 1 conditions; total gene numbers are shown in horizontal bar graphs, and comparison set numbers are displayed in vertical bar graphs. (G) Scatter plot comparing the log₂-fold changes for HS 1 (x-axis) and CS 1 (y-axis), color-coded according to plotting density; fold change densities for each condition are displayed in the side density plots. (H) TPM normalized expression values from the control condition for respective regulatory groups from panels A and D for CS 1 and HS 1 conditions; box edges delineate the first and third quartiles, center line represents the median, and whiskers denote points within 1.5× of the interquartile range.

Fig 3

Temporal analysis shows rapid rebalancing of transcriptome changes following heat shock, while cold shock elicits multiple responses. (A) Flow diagram visualizing the relative number and interconnected genes between conditions (x-axis) in the HS experiment for each regulatory group: strongly upregulated (dark green, bold font), upregulated (light green), non-regulated (white), strongly downregulated (dark brown, bold font), and downregulated (light brown). The flow diagram size is plotted to scale. Note that due to rounding of the percentages to the nearest whole number, some values might not add up to 100%. (B) Gene set enrichment analysis of archaeal clusters of orthologous genes (arCOGs) for all three HS conditions. Significance levels are indicated by a continuous color bar from white to dark purple. Genes with a P-value <0.05 are considered significantly overrepresented and highlighted by a circle, with circle size reflecting the number of genes in the category. For better comparison, up and down categories include all upregulated and downregulated genes, respectively, regardless of fold changes. (C) Flow diagram for the CS experiment displaying the relative number and interconnected genes between conditions (x-axis) for each regulatory group, as described in panel A. The flow diagram size is plotted to scale. Note that due to rounding of the percentages to the nearest whole number, some values might not add up to 100%. (D) Gene set enrichment analysis of arCOGs for all four CS conditions, with significance levels indicated by a continuous color bar from white to dark purple. Genes with a P-value <0.05 are considered significantly overrepresented and highlighted by a circle, with circle size reflecting the number of genes in the category. For better comparison, up and down categories include all upregulated and downregulated genes, respectively, regardless of fold changes.

Fig 4

Stress-induced genes are equipped with archaea-typical regulatory sequences features. (A) Position-specific promoter motifs of primary transcripts from (47), grouped by MEME motif search into classical archaeal promoter motif or not, based on the presence of BRE and TATA elements. (B) Transcript abundance (TPM-normalized counts under control conditions) for genes with a promoter motif (blue) and without (gray). Box edges delineate first and third quartiles, center line represents median, and whiskers indicate points within 1.5× interquartile range. Welch’s t-test used to assess differences between groups; significance level **** indicates P-values <0.0001. (C) Promoter strength comparison, estimated by P-values of MEME-detected promoter motifs (mode: one occurrence per site), for highly upregulated (dark green) and unchanged (white) genes under CS 1 and HS 1. Box plot parameters as in B; ns signifies P-values ≥0.05. (D) 5′ UTR length comparison for highly upregulated or downregulated genes against total distribution (Ctrl, 779). Density plots shown with bars in a window size of 1 and overlaid densities. (E) Proportion of leaderless (5′ UTR = 0) and leadered (5′ UTR ≥1) transcripts in highly upregulated and downregulated groups across all conditions. Chi-square test compares each stress condition to Ctrl sample distribution. Relative proportion color-coded from white (0) to dark purple (100%); significance (P < 0.05) indicated by an asterisk. (F) Position-specific terminator motif based on primary 3′ ends derived from Term-seq experiments, grouped by MEME search for genes with poly(U) signal (blue) or not (white). (G) Transcript abundance (TPM-normalized counts under control conditions) for genes with terminator poly(U) motif (blue) and without (gray). Box plot parameters as in B; significance level **** indicates P-values <0.0001. (H) Nucleotide enrichment meta-analysis compares nucleotide content of each group in a position-specific manner to randomly selected intergenic positions (n = 100,000). Log₂ enrichment shown for all primary 3′ ends with poly(U) motif (blue) or not (gray), and highly upregulated and downregulated genes in HS 1 and CS 1. (I) Proportion comparison of genes in regulatory groups with poly(U) signal or not, for all conditions. Note that an insufficient number of genes was detected for the CS 2 sample for robust statistical testing. Chi-square test checks whether the proportion of genes with poly(U) terminator differs significantly from Ctrl condition; no significant difference at P <0.05. (J) 3′ UTR comparison of regulatory groups to the Ctrl condition. Box plot parameters as in B; significance levels of * and ns indicate P-values <0.05 and ≥0.05, respectively.

Fig 5

Long-read sequencing reveals operon organization of heat shock genes. (A) Operon analysis using long-read PCR-cDNA Nanopore sequencing for HSP20 operon and (B) thermosome operon. Top panel: annotated operons from DOOR2 database, previous short-read RNA-sequencing, and long-read Ctrl condition sequencing. Genes visualized as rectangles, strand indicated by arrow direction, and operons connected by blue background lines. Next panel: current gene annotation and log₂-fold changes on RNA level for three HS conditions (1, 2, R) from left to right. Significance indicated by up or down arrows. Profiles of mean normalized coverage values from Nanopore and short-read RNA-sequencing, color-coded as Ctrl (gray) and HS 1 (orange). The last panel presents a single-read analysis of Nanopore reads. Each line represents one full-length sequenced read, sorted by their read start position. Vertical lines indicate primary start sites identified in reference (47), while dashed lines represent primary 3′ ends from short-read Term-sequencing.

Fig 6

Proteomics count values but not fold changes are moderatly correlated to RNA levels. (A) Volcano plot displaying log₂-fold changes in protein levels (x-axis) and −log10 (adjusted P-value, padj) on the y-axis comparing HS 1 (5 min) with control condition. Protein-coding genes are categorized based on significance and fold changes, color-coded as upregulated (padj < 0.05, log₂FC > 0, green), non-regulated (NS, padj ≥ 0.05, white), and downregulated (padj < 0.05 and log₂FC < 0, brown). HS genes are highlighted. Genes were assessed at a significance level of 0.05, indicated by a dashed line in the plot. (B) HS and (C) CS flow diagram visualizing the relative number and interconnected genes between conditions on the x-axis for each regulatory group, as described in panel A. Note that due to rounding of the percentages to the nearest whole number, some values might not add up to 100%. (D) Gene set enrichment analysis of archaeal clusters of orthologous genes (arCOGs) for all three HS and (E) all four CS conditions. Significance level is indicated by a continuous color bar from white to dark purple. Genes with a P-value <0.05 are considered significantly overrepresented in the category and are highlighted by a circle. The circle size reflects the number of genes in the category. (F) Comparison of log₂-fold changes and (G) normalized TPM and iBAQ values of HS 1 measured at RNA (x-axis) and protein (y-axis) levels. HS genes are highlighted. Pearson’s correlation coefficient is shown in the top left. Density of plotting is color-coded. (H) Correlation matrix (Pearson’s correlation coefficient) of pairwise comparisons based on fold changes (upper left corner) and normalized count values (bottom right).

Fig 7

Cluster analysis identifies signature genes involved in heat shock response. (A) Pathway plot illustrating log₂-fold changes at RNA level (x-axis) and protein level (y-axis) in a line plot from condition HS 1 via HS 2 to recovery condition. Genes of the Phr regulon are highlighted with colors and points. (B) Fold changes for genes presented in panel A. Significant regulation in respective conditions is indicated by a circle (padj < 0.05) or a rectangle. (C) Comparison of median z-score values shown as a point in each cluster (row) and condition (column) for protein (yellow) and RNA (gray) values. Shaded area represents the interquartile range. (D) Log₂-fold changes of genes sorted by clusters. Points indicate median, while shaded area shows interquartile range. (E) Gene set enrichment analysis of archaeal clusters of orthologous genes (arCOGs) for five selected clusters. Significance level is represented by a continuous color bar from white to dark purple. Genes with a P-value <0.05 are considered significantly overrepresented in the category and are highlighted by a circle. Circle size reflects the number of genes in the category. (F) Functional enrichment analysis of genes from each cluster was performed by selecting all significantly regulated genes in the most significantly overrepresented arCOG category. Confidence in predicted interactions, according to the STRING database, is indicated by line thickness. Only genes with at least one connection are shown. Genes are depicted as points and colored based on the fold change at RNA level in condition HS 1.

Fig 8

Cluster analysis identifies signature genes involved in cold shock response. (A) Pathway plot illustrating log₂-fold changes at RNA level (x-axis) and protein level (y-axis) in a line plot from condition CS 1, through CS 2 and CS 3 to recovery condition. The nine genes with the highest upregulation at RNA level (CS 1) are highlighted with colors and points. (B) Fold changes for genes presented in panel G. Significant regulation in respective conditions is indicated by a circle (padj < 0.05) or a rectangle. (C), Comparison of median Z-score values shown as a point in each cluster (row) and condition (column) for protein (yellow) and RNA (gray) values. Shaded area represents the interquartile range. (D) Log₂-fold changes of genes sorted by clusters. Points indicate median, while shaded area shows interquartile range. (E) Gene set enrichment analysis of archaeal clusters of orthologous genes (arCOGs) for four selected clusters. Significance level is represented by a continuous color bar from white to dark purple. Genes with a P-value <0.05 are considered significantly overrepresented in the category and are highlighted by a circle. Circle size reflects the number of genes in the category. (F) Functional enrichment analysis of genes from each cluster was performed by selecting all significantly regulated genes in the most significantly overrepresented arCOG category. Confidence in predicted interactions, according to the STRING database, is indicated by line thickness. Only genes with at least one connection are shown. Genes are depicted as points and colored based on the fold change at protein level in condition CS 3.