Synthetic messenger RNA vaccines and transcriptomic dysregulation: Evidence from new-onset adverse events and cancers post-vaccination

Abstract

BACKGROUND

Synthetic messenger RNA (mRNA) vaccines have raised concerns regarding prolonged spike protein expression, immune activation, and potential off-target effects.

AIM

To investigate transcriptomic alterations in individuals with new-onset adverse events or cancer following mRNA coronavirus disease 2019 vaccination.

METHODS

Bulk RNA sequencing was performed on peripheral blood from two patient groups: (1) Individuals with new-onset nonmalignant adverse events; and (2) Individuals newly diagnosed with cancer post-vaccination. A control group of normal individuals was used for comparison. Differential gene expression was analyzed using DESeq2, and Gene Set Enrichment Analysis was conducted using the MSigDB database and custom gene sets.

RESULTS

Both vaccine patient groups displayed widespread transcriptional dysregulation. In the nonmalignant adverse event group, hallmark enrichments included mitochondrial dysfunction, proteasome-mediated stress, transcriptomic instability, and systemic inflammation. The cancer group exhibited additional hallmarks of genomic instability and epigenetic reprogramming. Nonsense-mediated decay, ribosomal stress, and myelocytomatosis oncogene activation were prominent in both groups, while immune signaling via toll-like receptors and type I interferons was particularly elevated in cancer patients. The observed transcriptomic profiles indicate cellular stress responses, mitochondrial dysfunction, and immune dysregulation following exposure to mRNA vaccines, potentially in susceptible individuals.

CONCLUSION

Shared and distinct molecular signatures in both cohorts demonstrate underlying mechanisms contributing to post-vaccine symptomatology and complications, including oncogenesis and or progression of malignant disease. These findings underscore the need for a deeper investigation into the long-term safety of mRNA vaccines and host response variability.

Key Words: Coronavirus disease; Vaccine; RNA-seq; Immune dysregulation; Cancer

Core Tip: This study demonstrates that individuals experiencing new-onset adverse events or cancer after messenger RNA (mRNA) coronavirus disease 2019 vaccination exhibit widespread transcriptomic dysregulation. Bulk RNA sequencing revealed hallmarks of mitochondrial dysfunction, systemic inflammation, proteasome and ribosomal stress, and nonsense-mediated decay, with additional genomic instability and epigenetic reprogramming in cancer patients. Notably, myelocytomatosis oncogene activation and heightened immune signaling via toll-like receptors and type I interferons were observed. These findings highlight shared and distinct molecular signatures, underscoring the need for further investigation into long-term mRNA vaccine safety and host variability.

INTRODUCTION

Synthetic messenger RNA (mRNA) vaccines represent a novel immunization platform that delivers chemically modified mRNA, encapsulated in lipid nanoparticles (LNPs), to encode and express viral antigens in human cells. The accelerated development of mRNA vaccines was primarily driven by the urgency of the coronavirus disease 2019 (COVID-19) pandemic, supported by emergency authorizations, platform-based manufacturing, and the parallel rather than sequential conduct of clinical trial phases. While initially promoted for their rapid production and strong immunogenicity, a growing body of evidence has identified a range of adverse events, including myocarditis, thrombosis, cerebrovascular accidents, amyloidogenesis, arthralgia, menstrual disorders, reproductive concerns, and new-onset malignancies[1-5]. These events have prompted renewed scrutiny of the molecular mechanisms triggered by synthetic mRNA constructs and their intracellular fate[6,7].

Unlike endogenous mRNA, vaccine-derived transcripts incorporate non-natural features such as N1-methylpseudouridine (m1Ψ) substitution, extended poly(A) tails, and optimized untranslated regions to enhance stability and translational efficiency. These same modifications may alter RNA metabolism and surveillance, contributing to ribosomal infidelity, +1 frameshifting, and defective protein folding[8,9]. Independent studies have reported reverse transcription of vaccine mRNA into DNA via endogenous long interspersed nuclear element-1 (LINE-1) elements, highlighting issues over genomic integration and persistent expression[8,10]. Recent analyses of BNT162b2 vaccine vials revealed residual plasmid DNA, including SV40 enhancer and antibiotic resistance genes, raising concerns over potential oncogenic risk from unintended genomic exposure[11-15]. Moreover, recent research has shown that spike protein translated from mRNA vaccines can persist in the circulation for extended periods, far beyond initial pharmacokinetic expectations, potentially sustaining inflammatory signaling and immune activation[16-18]. Finally, some studies suggest that severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) bears molecular signatures consistent with synthetic genome assembly, raising further concerns about engineered RNA platforms and their potential unforeseen biological consequences[19]. These features collectively raise concerns regarding unintended cellular consequences of synthetic mRNA exposure.

To address this gap, we conducted a comparative transcriptomic analysis using bulk RNA-sequencing data from whole blood of patients with new-onset adverse events (group 1) and patients with new-onset cancers (group 2) following COVID-19 mRNA vaccination, alongside normal controls. Differential expression analysis and Gene Set Enrichment Analysis (GSEA) were performed to identify hallmark-specific molecular pathways perturbed in each group. By comparing these two profiles, this study aims to highlight potential transcriptomic disturbances associated with vaccine-related pathology and explore whether distinct or shared molecular hallmarks may underlie different post-vaccination clinical trajectories.

MATERIALS AND METHODS

Patient population characteristics and sample collection

We enrolled adult participants (≥ 18 years) who developed new-onset symptoms following administration of mRNA-based COVID-19 vaccines [BNT162b2 (pfizer) or mRNA-1273 (moderna)]. All participants provided written informed consent prior to inclusion.

The study population was divided into two groups: (1) Group 1: Consisting of three individuals who experienced post-vaccination non-malignant adverse events with symptoms emerging within one month of vaccination; and (2) Group 2: Comprising seven individuals diagnosed with new-onset malignancies within one year of vaccination. None of the participants in either group had a prior history of chronic disease or known genetic predisposition before vaccination and symptom onset. The clinical manifestations of each patient are summarized in Table 1. As a reference, RNA-seq data from 803 unvaccinated individuals were retrieved from the GTEx dataset[20]. These control samples were collected prior to the development of COVID-19 vaccines, ensuring that all control individuals were unvaccinated[21].

Table 1 Summary of study participants, clinical manifestations, and diagnostic groupings.

Patient ID	Group	Condition/symptoms
P1	1: Non-malignant adverse events	Brain fog and chronic fatigue
P2	1: Non-malignant adverse events	Cardiovascular injury, thrombosis, and chronic fatigue
P3	1: Non-malignant adverse events	Cardiovascular injury, thrombosis, and chronic fatigue
P4	2: Cancer	Glioblastoma multiforme
P5	2: Cancer	Prostate adenocarcinoma
P6	2: Cancer	Bladder tumor
P7	2: Cancer	Follicular lymphoma
P8	2: Cancer	Prostate cancer
P9	2: Cancer	Glioblastoma multiforme
P10	2: Cancer	Follicular lymphoma

Peripheral blood samples were collected by licensed nurses at Neo7Bioscience-affiliated clinical sites, following institutional biosafety protocols. Venipuncture was performed using standard sterile technique, and whole blood was drawn into pre-labelled streck tubes and stored at 4 °C. Samples from group 1 (post-vaccination symptom group) were transported the same day in cooled containers to the Genomics Center at the University of North Texas for processing. Samples from group 2 (post-vaccination malignancy group) were collected into PAXgene Blood RNA Tubes (Qiagen) and transported at controlled temperature to Psomagen-Multiomics services and data analysis for downstream processing.

RNA extraction

Total RNA was extracted from whole blood using column-based purification workflows compatible with downstream RNA-seq. For group 1, RNA was extracted with the Quick-DNA/RNA Viral Kit (Zymo Research), including on-column DNase treatment to remove genomic DNA contamination. For group 2, RNA was extracted from the entire blood volume of a single PAXgene Blood RNA Tube using the PAXgene Blood RNA Kit (Qiagen), following the manufacturer’s protocol. Extracted RNA was quantified and stored at -20 °C until further use. For both groups, RNA quality was assessed by electrophoretic profiling (Agilent TapeStation 4200 or equivalent), and only samples meeting quality thresholds (e.g., RNA integrity number ≥ 7.0 or equivalent QC metrics) were advanced to library preparation.

Library preparation and complementary DNA synthesis and high-throughput sequencing

RNA-seq libraries were generated using strand-specific protocols with ribosomal RNA (rRNA) depletion to maximize detection of coding transcripts. For group 1, libraries were prepared using the Illumina TruSeq Stranded Total RNA Kit with Ribo-Zero depletion chemistry (Illumina). RNA was enzymatically fragmented, reverse-transcribed to complementary DNA (cDNA), end-repaired, A-tailed, adapter-ligated, and polymerase chain reaction (PCR)-amplified. Libraries were quantified fluorometrically, fragment sizes were verified by electrophoresis, and sequencing was performed on an Illumina NextSeq 550 (High Output v2.5, 300 cycles) with paired-end 150 bp reads.

For group 2, libraries were prepared using the Illumina Stranded Total RNA Ribo-Zero Plus Kit (Illumina), incorporating rRNA depletion and strand-specific cDNA synthesis. RNA was enzymatically fragmented, reverse-transcribed to first-strand cDNA, converted to double-stranded cDNA, end-repaired, adenylated, ligated to Illumina-specific adapters, and PCR-amplified. Library quality and fragment size distribution were assessed using D5000 ScreenTape on the Agilent TapeStation 4200. Libraries were stored at -20 °C before sequencing. Normalized and pooled libraries were sequenced on the Illumina NovaSeq X platform (NovaSeq X Series 25B Reagent Kit, 300 cycles), generating 151 bp paired-end reads to a target yield of 3.02 Gb per sample (approximately 20 million reads). PhiX was added at 1% as an internal control, and sequencing quality metrics (Q30 scores, %PF) were monitored in real time using Illumina’s Sequencing Analysis Viewer.

Data preprocessing, quality control, and alignment

The raw sequencing data, in the form of paired-end fastq files, were first organized by sample and read type (R1 and R2). These fastq files were stored in a dedicated directory, and the cat command was used to group them based on their sample identifiers. The R1 and R2 files for each sample were then merged into a single file for each read type, resulting in forward and reverse merged files that were ready for downstream analysis.

Quality control checks were performed before the alignment step to ensure the integrity of the data. The sequencing reads were aligned to the human reference genome (hg38) using STAR[22], a widely used aligner for RNA sequencing data. STAR efficiently handles read alignment in a multi-threaded environment, optimizing processing time. The aligned reads were then outputted as BAM files, sorted by coordinate, which is standard practice for subsequent analysis steps.

Differential expression analysis and volcano plot visualization

To investigate transcriptomic alterations associated with mRNA vaccine–related outcomes, we performed two distinct differential expression analyses, each comparing one patient group to a common pool of normal controls (n = 803)[20]. The first comparison included individuals from group 1, who developed nonmalignant new-onset symptoms following mRNA vaccination (n = 3), while the second comprised individuals from group 2, diagnosed with cancer as a new-onset condition shortly after vaccination (n = 7). For both analyses, raw gene counts generated from Salmon[23] were input into DESeq2[24], which applied internal normalization and modeled dispersion using the negative binomial distribution. Shrinkage of log₂ fold changes was performed using the “apeglm” method to improve effect size estimation, particularly for genes with low read counts.

To visualize the global patterns of gene expression, volcano plots were generated for each comparison. In these plots, the X-axis represents the log₂ fold change, and the Y-axis represents the negative log₁₀ of the adjusted P value. Genes with log₂ fold change greater than +1 and an adjusted P value (P_adj) below 0.05 were considered significantly upregulated and were colored red, while those with log₂ fold change less than -1 and P_adj < 0.05 were considered significantly downregulated and were colored blue. All remaining genes that did not meet these criteria were displayed in gray.

GSEA

To investigate transcriptomic disruptions associated with mRNA vaccine exposure, we analyzed RNA-seq data from two case groups: (1) Group 1: Comprising three individuals who developed new-onset vaccine-related adverse effects; and (2) Group 2: Consisting of seven individuals diagnosed with new-onset cancer following mRNA COVID-19 vaccination. Both groups were compared to a shared reference cohort of 803 individuals in a normal control group, whose RNA-seq data were retrieved from the GTEx dataset[20]. Differential gene expression analysis was conducted using DESeq2[24], and genes were ranked by log₂ fold change to generate the input file for enrichment analysis.

We then performed GSEA using the preranked mode of the Broad Institute’s GSEA 4.4.0 tool[25]. This approach was selected due to the significant class imbalance (3 samples vs 803 samples), as it improves the stability and interpretability of enrichment results. The “weighted” enrichment statistic was applied to incorporate both gene ranking and expression magnitude, enhancing sensitivity to biologically relevant perturbations. We used 1000 gene set permutations to estimate statistical significance and selected the “no collapse” option to retain gene symbols without alias mapping.

GSEA was systematically conducted across multiple MSigDB[26] collections: (1) H (Hallmark); (2) C2 (Kyoto Encyclopedia of Genes and Genomes and REACTOME); (3) C5 (Gene Ontology Biological Processes); (4) C7 (immunologic signatures); and (5) Custom gene sets including Gavish and curated grape seed extract-derived expression signatures. Enrichment results were interpreted based on the Normalized Enrichment Score (NES) and false discovery rate (FDR) (q value), selecting gene sets that were strongly enriched (positive NES) or suppressed (negative NES), which indicated upregulated or downregulated biological functions, respectively. For downstream analysis, only gene sets with FDR < 0.25 and absolute NES ≥ 1.5 were retained. The top enriched and suppressed gene sets were manually grouped into higher-order molecular pathway categories related to the patient’s group symptoms.

To improve biological relevance, we excluded gene sets associated with non-blood tissues (e.g., brain, retina, skin) and embryonic development, as these are unlikely to reflect transcriptional activity in peripheral blood. This filtering minimized noise and focused the analysis on pathways relevant to immune, inflammatory, and systemic responses.

Protein-protein interaction map

For constructing protein-protein interaction (PPI) networks, interaction information of differentially expressed genes from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database was loaded onto the Cytoscape software platform[27,28]. The STRING database is a comprehensive online resource that compiles and predicts PPIs across numerous organisms[29,30]. It integrates known and predicted associations derived from experimental data, computational prediction methods, co-expression analyses, text mining, and curated databases. Each interaction is scored and visualized in an interactive network, helping researchers explore functional associations between proteins within cellular pathways or biological processes[31]. The network can be viewed within STRING or can be exported directly to Cytoscape. Cytoscape combines biomolecular interaction networks with high-throughput expression data and other molecular states into an integrated conceptual network model.

RESULTS

To characterize the global transcriptional impact of mRNA vaccination in individuals with post-vaccination adverse outcomes, we conducted differential gene expression analyses comparing each patient group to the normal control cohort. As shown in Figure 1, the volcano plots visualize the distribution of gene expression changes for both group 1, new-onset nonmalignant post-vaccination symptoms (Figure 1A), and group 2, new-onset cancer diagnoses (Figure 1B). In both groups, a clear transcriptomic shift is observed, with hundreds of genes showing significant dysregulation. Notably, the cancer group exhibits a broader distribution of downregulated genes, including several with extreme fold change magnitudes, suggesting a profound suppression of key regulatory pathways.

Figure 1 Transcriptome-wide differential expression profiles in vaccine-affected individuals. A: Volcano plot showing differential gene expression in individuals with new-onset adverse events following messenger RNA (mRNA) coronavirus disease 2019 (COVID-19) vaccination (n = 3) compared to normal controls (n = 803); B: Volcano plot showing differential gene expression in individuals diagnosed with new-onset cancer shortly after receiving mRNA COVID-19 vaccination (n = 7) compared to the same control cohort. Each point represents a single gene plotted by log₂ fold change (X-axis) and -log₁₀ adjusted P value (P_adj) (Y-axis). Genes with significant upregulation (log₂FC > 1, P_adj < 0.05) are marked in red, while significantly downregulated genes (log₂FC < -1, P_adj < 0.05) are shown in blue. Non-significant genes appear in gray. These plots reveal widespread transcriptional dysregulation in both patient groups, serving as the foundation for subsequent pathway enrichment analysis.

GSEA was performed separately for both study groups: (1) Group 1: Comprising individuals with new-onset adverse events following mRNA COVID-19 vaccination; and (2) Group 2: Comprising individuals diagnosed with new-onset cancer after vaccination, to characterize transcriptomic disruptions associated with vaccine exposure. GSEA was conducted using the predefined MSigDB gene set collections, and the enriched gene sets were systematically grouped according to their relevance to higher-order molecular hallmark pathways (Tables 2 and 3). This pathway-based grouping enabled integrative interpretation of transcriptomic alterations across systems biology domains.

Table 2 Gene set enrichment results for group 1 categorized by molecular pathway.

Molecular pathways	MSigDB standard name	Normalized Enrichment Score	False discovery rate
Mitochondrial electron transport dysfunction and reactive oxygen species	KEGG_MEDICUS_REFERENCE_MITOCHONDRIAL_COMPLEX_UCP1_IN_THERMOGENESIS	1.86	0.11
	KEGG_MEDICUS_REFERENCE_ELECTRON_TRANSFER_IN_COMPLEX_I	1.70	0.10
	KEGG_MEDICUS_VARIANT_MUTATION_INACTIVATED_PINK1_TO_ELECTRON_TRANSFER_IN_COMPLEX_I	1.73	0.09
	KEGG_MEDICUS_VARIANT_MUTATION_CAUSED_ABERRANT_ABETA_TO_ELECTRON_TRANSFER_IN_COMPLEX_I	1.69	0.09
	KEGG_MEDICUS_VARIANT_MUTATION_CAUSED_ABERRANT_SNCA_TO_ELECTRON_TRANSFER_IN_COMPLEX_I	1.66	0.11
	HALLMARK_OXIDATIVE_PHOSPHORYLATION	1.43	> 0.001
Proteasome-mediated protein degradation stress	KEGG_MEDICUS_VARIANT_MUTATION_CAUSED_ABERRANT_ABETA_TO_26S_PROTEASOME_MEDIATED_PROTEIN_DEGRADATION	1.78	0.12
	KEGG_MEDICUS_VARIANT_MUTATION_INACTIVATED_VCP_TO_26S_PROTEASOME_MEDIATED_PROTEIN_DEGRADATION	1.74	0.13
	KEGG_MEDICUS_VARIANT_MUTATION_CAUSED_ABERRANT_SOD1_TO_26S_PROTEASOME_MEDIATED_PROTEIN_DEGRADATION	1.73	0.09
	GOBP_PROTEIN_CONTAINING_COMPLEX_REMODELING	-2.11	> 0.001
Transcriptomic instability and translational stress	REACTOME_EUKARYOTIC_TRANSLATION_INITIATION	3.11	> 0.001
	REACTOME_SRP_DEPENDENT_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE	3.08	> 0.001
	REACTOME_NONSENSE_MEDIATED_DECAY_NMD	2.9	> 0.001
	KEGG_MEDICUS_REFERENCE_TRANSLATION_INITIATION	3.0	> 0.001
	GOBP_CYTOPLASMIC_TRANSLATION	2.60	> 0.001
	GOBP_RIBOSOMAL_SMALL_SUBUNIT_BIOGENESIS	2.35	> 0.001
	GOBP_RIBOSOMAL_LARGE_SUBUNIT_BIOGENESIS	2.34	> 0.001
Systemic inflammatory and immune response	GSE22886_NAIVE_TCELL_VS_DC_UP	2.43	> 0.001
	GSE2405_0H_VS_9H_A_PHAGOCYTOPHILUM_STIM_NEUTROPHILS_UP	2.36	> 0.001
	GSE7400_CTRL_VS_CSF3_IN_VIVO_TREATED_PBMC_UP	2.27	> 0.001
	GSE7400_CTRL_VS_CSF3_IN_VIVO_TREATED_PBMC_DN	-1.93	> 0.001
	GSE6269_E_COLI_VS_STREP_PNEUMO_INF_PBMC_DN	-1.9	> 0.001
	GOBP_ANTIMICROBIAL_HUMORAL_RESPONSE	-1.95	> 0.001
	GOBP_HUMORAL_IMMUNE_RESPONSE	-1.91	0.006
Endothelium dysfunction	HALLMARK_ANGIOGENESIS	-1.94	> 0.001
Proliferative signaling and suppressed tumor control	HALLMARK_MYC_TARGETS_V1	1.64	> 0.001
	GAVISH_3CA_MALIGNANT_METAPROGRAM_3_CELL_CYLCE_HMG_RICH	1.77	0.005
	HALLMARK_KRAS_SIGNALING_DN	-1.68	> 0.001
	HALLMARK_P53_PATHWAY	-1.50	0.038
	KEGG_MEDICUS_REFERENCE_WNT_SIGNALING_MODULATION_WNT_INHIBITOR	-1.56	0.89

The Table 2 displays the top gene sets within each category, along with their corresponding MSigDB identifiers, Normalized Enrichment Score (NES), and false discovery rate (FDR) q value. Positive NES values indicate upregulation, while negative NES values indicate downregulation of the gene sets in vaccinated samples (group 1) relative to controls. Gene sets were selected based on biological relevance and NES >|15| with FDR q < 0.25. Group 1 includes patients with new-onset adverse events following messenger RNA coronavirus disease 2019 vaccination (n = 3), compared with normal controls (n = 803 unvaccinated individuals from the GTEx dataset).

Table 3 Gene set enrichment results for group 2 categorized by molecular pathway.

Molecular pathways	MSigDB standard name	Normalized Enrichment Score	False discovery rate
Transcriptomic instability and translational stress	REACTOME_SIRT1_NEGATIVELY_REGULATES_RRNA_EXPRESSION	3.08	> 0.001
	REACTOME_RNA_POLYMERASE_I_PROMOTER_ESCAPE	2.61	> 0.001
	REACTOME_TRANSCRIPTIONAL_REGULATION_BY_SMALL_RNAS	2.63	> 0.001
	REACTOME_B_WICH_COMPLEX_POSITIVELY_REGULATES_RRNA_EXPRESSION	2.60	> 0.001
	REACTOME_POSITIVE_EPIGENETIC_REGULATION_OF_RRNA_EXPRESSION	2.54	> 0.001
	REACTOME_RRNA_MODIFICATION_IN_THE_NUCLEUS_AND_CYTOSOL	2.03	0.027
Systemic inflammatory and immune response	KEGG_MEDICUS_REFERENCE_TYPE_I_INTERFERON_TO_JAK_STAT_SIGNALING_PATHWAY	2.64	> 0.001
	KEGG_MEDICUS_PATHOGEN_HIV_TAT_TO_TLR2_4_NFKB_SIGNALING_PATHWAY	2.73	> 0.001
	KEGG_MEDICUS_REFERENCE_TYPE_I_IFN_SIGNALING_PATHWAY	2.64	> 0.001
	KEGG_MEDICUS_REFERENCE_RIG_I_NFKB_SIGNALING_PATHWAY	2.54	> 0.001
	KEGG_MEDICUS_REFERENCE_TLR3_IRF7_SIGNALING_PATHWAY	2.46	> 0.001
	KEGG_MEDICUS_REFERENCE_TLR7_8_9_IRF5_SIGNALING_PATHWAY	2.29	> 0.001
	KEGG_MEDICUS_REFERENCE_TLR7_9_IRF7_SIGNALING_PATHWAY	2.27	> 0.001
	KEGG_MEDICUS_REFERENCE_MDA5_IRF7_3_SIGNALING_PATHWAY	2.0	0.006
	KEGG_MEDICUS_REFERENCE_IFN_RIPK1_3_SIGNALING_PATHWAY	2.0	0.006
	KEGG_MEDICUS_REFERENCE_RIG_I_IRF7_3_SIGNALING_PATHWAY	2.0	0.007
	GAVISH_3CA_METAPROGRAM_CD8_T_CELLS_CHROMATIN	1.78	0.068
Endothelium dysfunction	GOBP_NEGATIVE_REGULATION_OF_ENDOTHELIAL_CELL_PROLIFERATION	-1.98	0.0017
	GOBP_NEGATIVE_REGULATION_OF_COAGULATION	-1.98	0.0017
Proliferative signaling and suppressed tumor control	REACTOME_ASSEMBLY_OF_THE_ORC_COMPLEX_AT_THE_ORIGIN_OF_REPLICATION	2.93	> 0.001
Genomic instability and epigenetic shift	REACTOME_DNA_METHYLATION	3.18	> 0.001
	REACTOME_CONDENSATION_OF_PROPHASE_CHROMOSOMES	2.93	> 0.001
	KEGG_MEDICUS_REFERENCE_CGAS_STING_SIGNALING_PATHWAY	2.43	> 0.001
	GAVISH_3CA_MALIGNANT_METAPROGRAM_4_CHROMATIN	1.91	0.05
	GOBP_NUCLEOSOME_ORGANIZATION	2.33	0.01

The Table 3 displays the top gene sets within each category, along with their corresponding MSigDB identifiers, Normalized Enrichment Score (NES), and false discovery rate (FDR) q value. Positive NES values indicate upregulation, while negative NES values indicate downregulation of the gene sets in vaccinated samples relative to controls. Gene sets were selected based on biological relevance and NES >|15| with FDR q < 0.25. Group 2 includes patients with new-onset cancers diagnosed shortly after messenger RNA coronavirus disease 2019 vaccination (n = 7), compared with normal controls (n = 803 unvaccinated individuals from the GTEx dataset).

In group 1, the top enriched gene sets clustered into six key hallmark categories: (1) Mitochondrial electron transport dysfunction and reactive oxygen species (ROS); (2) Proteasome-mediated protein degradation stress; (3) Transcriptomic instability and translational stress; (4) Systemic inflammatory and immune response; (5) Endothelium dysfunction; and (6) Proliferative signaling and suppressed tumor control. In contrast, group 2 exhibited a distinct but overlapping enrichment pattern. The hallmarks identified in this group include: (1) Transcriptomic instability and translational stress; (2) Systemic inflammatory and immune response; (3) Endothelium dysfunction; (4) Proliferative signaling and suppressed tumor control; and (5) Genomic instability and epigenetic shift. As shown in Table 1, multiple gene sets within each hallmark category demonstrated strong enrichment scores (NES > 1.5 or < -1.5) and statistically significant FDR q values (FDR < 0.25), supporting the presence of coordinated molecular dysregulation. Representative enrichment score curves are provided in Supplementary Figures 1 and 2 to illustrate the distribution of gene ranks contributing to each hallmark. In addition, Supplementary Tables 1 and 2 provide the systematic names, standard names, and external links or source publications for all enriched gene sets.

To further explore the molecular landscape associated with vaccine-induced transcriptomic disruptions, PPI network analysis was conducted for the most significantly dysregulated genes in each study group. As shown in Figure 2, the PPI map for group 1 (individuals with new-onset adverse events post-vaccination) and Figure 3 shows the PPI map for group 2 (individuals with new-onset cancer following vaccination)

Figure 2 Protein-protein interaction network of the most dysregulated genes in group 1. A: Genes related to proliferative signaling and tumor control; B: Genes related to systemic inflammatory immune; C: Genes related to transcriptomic instability; D: Genes related to mitochondrial electron transport dysfunction and reactive oxygen species; E: Genes related to proteasome-mediated protein degradation stress. Group 1 includes patients with new-onset adverse events following messenger RNA coronavirus disease 2019 vaccination (n = 3), compared with normal controls (n = 803 unvaccinated individuals from the GTEx dataset). Node color intensity reflects the degree of interaction (connectivity), with darker nodes indicating higher connectivity or a hub status within the network.

Figure 3 Protein-protein interaction network of the most dysregulated genes in group 2. A: Genes related to transcriptomic instability, translational stress; B: Genes related to genomic instability and epigenetic shift; C: Genes associated with endothelial dysfunction; D: Genes related to systemic inflammatory and immune response; E: Genes that are related to proliferative signaling and suppressed tumor control. Group 2 includes patients with new-onset cancers diagnosed shortly after messenger RNA coronavirus disease 2019 vaccination (n = 7), compared with normal controls (n = 803 unvaccinated individuals from the GTEx dataset). Node color intensity reflects the degree of interaction (connectivity), with darker nodes indicating higher connectivity or a hub status within the network.

DISCUSSION

The rapid global rollout of RNA-based COVID-19 vaccines introduced a novel therapeutic platform involving synthetic mRNA and LNP delivery systems. While initially deployed to mitigate the spread of SARS-CoV-2, the long-term biological effects were unknown. Unlike conventional vaccines, these formulations induce host cells to express a viral spike glycoprotein from exogenous mRNA templates, raising concerns about unintended cellular responses[32]. Emerging clinical reports of persistent symptoms following vaccination – including neurological, cardiovascular, and immunological disturbances – have prompted scientific scrutiny into the molecular mechanisms potentially triggered by prolonged exposure to synthetic mRNA, its translation products, and associated immune activation[2,3,33,34]. To investigate these concerns, we performed GSEA on transcriptomic data from two distinct patient cohorts: (1) Individuals with new-onset post-vaccination adverse effects (group 1); and (2) Patients who developed cancer following mRNA vaccination (group 2). A summary of the transcriptomic alterations observed in each group is illustrated in Figure 4.

Figure 4 Transcriptomic alterations following messenger RNA vaccination. This central illustration summarizes the experimental design, core findings, and proposed molecular mechanisms underlying transcriptomic dysregulation following synthetic messenger RNA (mRNA) coronavirus disease 2019 vaccination. Top left: Schematic of mRNA vaccination, showing lipid nanoparticle-encapsulated, chemically modified mRNA encoding spike protein delivered into host cells, initiating persistent translation and immunologic engagement. Top right: Volcano plots depict global differential gene expression in peripheral blood samples from two affected cohorts vs normal controls (n = 803). Left: Individuals with new-onset nonmalignant adverse events (n = 3). Right: Individuals with new-onset cancer (n = 7). Upregulated genes [log₂FC > 1, adjusted P value (P_adj) < 0.05] are shown in red; downregulated genes (log₂FC < -1, P_adj < 0.05) in blue; non-significant genes in gray. Bottom left (new-onset adverse events): Transcriptomic analysis reveals enrichment of pathways linked to mitochondrial electron transport dysfunction and reactive oxygen species, proteasome-mediated protein degradation stress, mRNA surveillance activation, and systemic inflammatory signaling. Bottom right (new-onset cancer): Cancer patients exhibit hallmarks of oncogenesis, including genomic instability, epigenetic reprogramming, nonsense-mediated decay, ribosomal stress, myelocytomatosis oncogene-driven proliferative signaling, and persistent immune activation via toll-like receptors and type I interferons. COVID-19: Coronavirus disease 2019; mRNA: Messenger RNA; MYC: Myelocytomatosis oncogene; ROS: Reactive oxygen species; TCA: Tricarboxylic acid; TLRs: Toll-like receptors;

In the group 1 cohort, enrichment of gene sets related to mitochondrial electron transport indicates transcriptional disruption in core components of oxidative phosphorylation, particularly at complex I. The presence of variant-related enrichments involving PTEN-induced putative kinase 1, amyloid-beta, and alpha-synuclein suggests convergence associated with impaired mitophagy, neurodegeneration, and defective electron flow. These findings suggest that spike-mediated inflammation disrupts mitochondrial homeostasis. Indeed, previous studies indicate that the spike protein alters mitochondrial fusion-fission dynamics, suppressing the biogenesis of its regulators (nuclear respiratory factor 1/2, peroxisome proliferator-activated receptor gamma coactivator 1-alpha, mitochondrial transcription factor A), and increasing mitochondrial damage[35]. In addition, numerous studies have independently reported persistent mitochondrial dysfunction as a hallmark of long COVID[36-38], further supporting a mechanistic link between spike protein exposure and bioenergetic imbalance. Mitochondrial electron transport and ROS-related shifts are depicted in the PPI network (Figure 2D), highlighting key interconnected hub genes.

In the group 1 cohort, enrichment of proteasome-associated gene sets indicated activation of components of the ubiquitin-proteasome system. Upregulated proteins such as ubiquitin-40S ribosomal protein S27a (RPS27A), 26S proteasome subunit SEM1, polyubiquitin-B (UBB), and ubiquitin-60S ribosomal protein L40 (UBA52), along with several core subunits of the 20S and 19S proteasome complexes, including PSMA2-5, PSMA8, PSMC1, and PSMD7 reflect increased transcriptional demand for protein degradation machinery, likely in compensation for misfolded or aggregated proteins. The enrichment of variant-associated gene sets linked to valosin-containing protein and super oxide dismutase 1 dysfunction further supports convergence with proteostasis disruption, a phenomenon commonly observed in neurodegenerative and proteinopathy-associated conditions. One potential driver of this sustained proteasomal activation is the prolonged presence of synthetic mRNA in circulation, which has been shown to persist beyond initial translation and may continuously stimulate the production of spike protein or aberrant translation products[16-18]. Moreover, persistence of vaccine-derived RNA has been demonstrated in human samples, including placental tissue up to 10 days post-vaccination[39] and plasma up to 28 days[40,41]. Identified persistence of vaccine-derived RNA in human heart tissue up to 30 days post-vaccination[42], demonstrated retention within lymph nodes for as long as 60 days, and Ota et al[43] recently reported mRNA detection in cerebral arteries 17 months post-vaccination[44]. Collectively, these findings underscore that synthetic mRNA and its byproducts may linger in diverse tissues for prolonged periods, continually engaging proteostatic and immune surveillance pathways. Additionally, emerging evidence suggests that the SARS-CoV-2 spike protein contains prion-like domains that are prone to misfolding and aggregation. These properties may be amplified in vaccine-induced expression, leading to persistent proteotoxic stress and neurodegenerative-like transcriptional profiles[8].

Transcriptomic Instability and Translational Stress was a dominant hallmark in both patient groups. mRNA stabilization, translation initiation, and protein synthesis rates are distinct yet interconnected regulatory nodes. The mRNA stabilization controls the substrate's half-life, determining the window of opportunity for translation. In contrast, translation initiation, often the rate-limiting step, involves the assembly of the ribosome and initiation factors (e.g., eukaryotic initiation factor 4F) at the 5' cap and is highly responsive to cellular signaling and stress. These processes can be uncoupled; a long-lived, stable mRNA may be poorly translated if initiation is blocked, while rapid synthesis can deplete short-lived mRNAs. Ultimately, the overall protein output is a product of both the mRNA's availability and the efficiency of the translational machinery[45-47]. Enrichment of gene sets related to translation initiation, ribosome biogenesis, and mRNA surveillance suggests sustained activation of the protein synthesis machinery. In group 1, upregulated pathways such as translation initiation, cytoplasmic translation, and ribosomal subunit biogenesis indicate persistent engagement of ribosomes and translational apparatus, likely driven by prolonged spike protein expression from stabilized synthetic mRNA[8,16,18]. Enrichment of nonsense-mediated decay (NMD) related genes in patients with new-onset vaccine-associated adverse events contrasts with previous findings that SARS-CoV-2 virus suppresses NMD to protect its RNA genome[48,49]. This suggests that, unlike the viral genome, synthetic mRNA used in vaccination may instead provoke a compensatory activation of RNA surveillance mechanisms, potentially due to persistent translation or accumulation of aberrant transcripts. The opposing patterns of NMD regulation in these two contexts point to distinct cellular responses and warrant further investigation. In group 2, as in the first group, we observed significant enrichment in gene sets related to ribosome biogenesis, cytoplasmic translation, and mRNA surveillance. However, the cancer group displayed more pronounced activation of ribosomal stress pathways, particularly those linked to RNA Polymerase I promoter escape, rRNA maturation, and positive epigenetic regulation of rRNA expression, indicating hyperactivation of nucleolar functions and elevated translational throughput. While these transcriptional programs are classically associated with tumor cells to support uncontrolled proliferation[50-52], their detection in peripheral blood samples likely reflects systemic consequences of underlying malignancy, such as systemic immune alterations or stress responses induced by tumor-related inflammation and signaling[53-55]. While the vaccine-adverse event group showed activation of NMD pathways, the cancer group displayed additional signals related to small RNA–mediated transcriptional regulation and epigenetic repression of ribosomal gene activity. The coexistence of transcriptional silencing and increased rRNA production may reflect underlying transcriptional stress affecting cellular balance.

Enrichments related to the hallmark of systemic inflammatory and immune response were identified in both groups analyzed in this work. In group 1, the most highly connected upregulated genes in the PPI network included LOX, CD28, CCR7, and SELL. Recent findings demonstrate that the use of m1Ψ in mRNA constructs can induce +1 ribosomal frameshifting, resulting in the production of off-target proteins that may elicit unintended cellular immune responses[9]. This aberrant antigen production may represent a novel mechanism contributing to systemic inflammation and immune dysregulation. In addition, emerging evidence suggests that cytoplasmic fragmentation of vaccine-derived mRNA may generate short RNA sequences with miRNA-like properties capable of hybridizing host immune transcripts, such as interferons and anti-inflammatory regulators. This unintended post-transcriptional interference could contribute to systemic immune dysregulation and inflammatory responses, particularly in individuals with predisposing comorbidities or impaired RNA degradation pathways[56]. In addition, numerous studies have previously demonstrated the inflammatory activity of mRNA LNPs[57-59]. Compared to group 1, where immune imbalance was more skewed toward aberrant cellular activation and humoral suppression, group 2 displayed a distinct pattern of innate immune activation. GSEA revealed the upregulation of pathways involving type I interferons, toll-like receptors, and nuclear factor kappa B (NF-κB)-driven inflammatory signaling. Enrichment of the RIG-I/MDA5–IRF7 axis, TLR3/7/8/9–IRF5/7, and the JAK-STAT interferon cascade suggests persistent engagement of RNA-sensing mechanisms and their downstream proinflammatory transcriptional programs. While these responses are central to antiviral immunity, their chronic activation in the peripheral blood of cancer patients might be related to inflammation, immune exhaustion, and tumor immune editing[60-62]. These immune alterations are consistent with previous findings in both elderly vaccine recipients and autoimmune patients, where IFN-JAK-STAT overactivation and RIG-I signaling were repeatedly identified as dominant signatures, underscoring that our observed systemic inflammation reflects a reproducible pattern across vaccinated cohorts[63].

For group 1, a negative NES for the angiogenesis hallmark suggests transcriptional downregulation or post-transcriptional inhibition of key angiogenic mediators. This may be linked to the RNA fragmentation hypothesis supported by Demongeot and Fougère’s work[56], wherein miRNA-like fragments derived from cleaved vaccine mRNA hybridize with endothelial transcripts, suppressing their translation and contributing to systemic endothelial dysfunction and impaired vascular repair. In the cancer group, transcriptomic analysis revealed significant downregulation of gene sets involved in negative regulation of endothelial cell proliferation and negative regulation of coagulation. The suppression of these regulatory pathways may indicate uncontrolled endothelial activation or a vascular pro-thrombotic shift, conditions that can be linked to spike protein adverse events and tumor progression[64-67]. Proteome-mediated protein degradation stress, marked by impaired ubiquitin-proteasome and autophagy pathways, contributes to the accumulation of misfolded or damaged proteins that sustain inflammation and cellular dysfunction. This ongoing proteotoxic stress is thought to contribute to the chronic fatigue, neurological dysfunction, and multi-organ symptoms in vaccinated patients[68,69].

For both groups 1 and 2, we could note an enrichment in the Proliferative Signaling and Suppressed Tumor Control Hallmark genes signatures. For group 1, positive enrichment in gene sets regulated by myelocytomatosis oncogene (MYC) suggests an active oncogenic transcriptional program that favors tumor growth and adaptation. This proliferative bias is compounded by the downregulation of tumor-suppressive pathways, including the KRAS-inhibited signature, p53 regulatory networks, and inhibitors of the Wnt pathway. Figure 2A depicts the PPI network for these hallmarks in group 1, highlighting the key interconnected genes driving these processes. These findings suggest a microenvironment conducive to unchecked cellular expansion, loss of apoptotic surveillance, and diminished responsiveness to anti-proliferative signals. Notably, many previous studies have reported that the spike protein can trigger the mitogen pathway, through the downregulation of angiotensin-converting enzyme 2 (ACE2) expression, which promotes an angiotensin II type I receptor (AT1R)-mediated signaling cascade, inducing the transcriptional regulatory molecules NF-κB and activator protein 1/c-Fos via mitogen-activated protein kinase activation[70-72]. In accordance with this molecular pathway, our investigated patients from group 1 and group 2 also presented a very downregulated ACE2 expression compared to normal control (logfold change: -4.3 and -4.8, respectively). In relation to this hallmark, group 2 showed enrichment of pathways involved in the DNA replication process, specifically the Assembly of the ORC Complex at the Origin of Replication, reflecting heightened replicative stress and uncontrolled proliferation potentially driven by tumor-related systemic effects or prolonged immune activation[60-62].

Only group 2 presented enrichments related to the Genomic Instability and Epigenetic Shift hallmark. Key pathways enriched show aberrant regulation of histone modification, DNA packaging, and epigenetic silencing, which are commonly observed in oncogenic events. In parallel, enrichment of the cGAS–STING signaling pathway points to innate immune recognition of cytoplasmic DNA fragments, a well-established marker of DNA damage and chromosomal instability[73]. This pathway is associated with tumor-promoting inflammation and immune editing[60-62]. Importantly, a recent longitudinal study of mRNA vaccination in octogenarians demonstrated that BNT162b2 administration elicited activation of the cGAS–STING pathway alongside robust immune and antibody responses, further supporting our transcriptomic findings[63]. Persistent genomic instability increases the likelihood of acquiring somatic mutations in key oncogenes and tumor suppressor genes, a process that can cumulatively drive malignant transformation. This is consistent with preclinical and in vitro evidence showing that exposure to the BNT162b2 mRNA vaccine can modulate endogenous reverse transcriptase activity (LINE-1), facilitate reverse transcription of vaccine mRNA into DNA, and alter nuclear localization of LINE-1 proteins[74], such events that could, in principle, contribute to insertional mutagenesis and genomic perturbations over time. Independent analyses have corroborated these concerns by identifying residual plasmid DNA contamination in both Pfizer and Moderna mRNA vaccines[75].

CONCLUSION

This study provides transcriptomic evidence of molecular disruptions in two patient populations – those with new-onset nonmalignant adverse events and those with newly diagnosed cancers. Using differential expression analysis and GSEA, we identified hallmark signatures of mitochondrial dysfunction, translational stress, immune dysregulation, endothelial disturbance, and proliferative signaling across both cohorts. Notably, while both groups shared transcriptional perturbations in immune and translational pathways, the cancer group exhibited additional signatures of genomic instability and epigenetic remodeling.

Persistent spike protein expression, prolonged synthetic mRNA activity, and RNA modifications such as m1Ψ appear to contribute to sustained aberrant ribosomal activity, proteostasis stress, and immune activation. Our findings also highlight transcriptional signals indicative of tumor-promoting conditions, including suppressed p53 networks, activated MYC targets, and altered interferon signaling, particularly in the context of epigenetic dysregulation in the cancer cohort. These observations suggest that vaccine-induced transcriptomic reprogramming may differentially affect individuals, genetically or immunologically, over a long period of time after vaccination. Despite the smaller sample size due to resource constraints (group 1 n = 3 vs n = 7 in group 2), the findings open an important avenue for understanding post-vaccine biological responses and underscore the value of expanding future studies with larger cohorts.

ACKNOWLEDGEMENTS

We thank the Genomics Center at the University of North Texas for sequencing support, Neo7Bioscience for data interpretation assistance, and the McCullough Foundation for scientific collaboration. We also acknowledge the participants who made this study possible.