Blog

Introduction

Rivaroxaban, the first orally active direct Factor Xa (FXa) inhibitor, is a small-molecule oxazolidinone derivative (a molecular weight of 435.88 Da). It specifically and reversibly binds to the S1 and S4 pockets of FXa with 10,000 times more selectivity than any other related serine proteases, such as thrombin, trypsin, plasmin, Factor VIIa, IXa, XIa, urokinase, or activated protein.^1,2 Rivaroxaban inhibits endogenous FXa activity with an inhibitory concentration 50% (IC₅₀) of 21±1 nM in a concentration-dependent manner¹ and thrombin generation is almost completely suppressed at therapeutic concentrations (80–100 nM).³

Based on previous landmark large Phase III randomized clinical trials (RCTs),^4–10 rivaroxaban is recommended for three indications: 1) stroke prevention of atrial fibrillation (AF), 2) reduction of major cardiovascular events in chronic coronary artery or peripheral artery disease, and 3) treatment and secondary prevention of venous thromboembolism (VTE) in adults and pediatric pupulations.^11–14

Since the approval of the US Food and Drug Administration (FDA) approval in 2011, the prescription volume of rivaroxaban has increased annually and reached 8.6 million in 2020, making it one of the most commonly prescribed medications in the United States.¹⁵ In China, rivaroxaban was approved for clinical use in 2009, and was the most frequently used direct oral anticoagulant (DOACs) in 2017.¹⁶

While demonstrating generally predictable pharmacokinetic (PK) and pharmacodynamic (PD) profiles, rivaroxaban exhibits substantial interindividual variability in therapeutic responses observed in both healthy populations and patient cohorts. The data from published studies have shown an up to 15-fold variation in rivaroxaban plasma concentration.¹⁷ During phase III trials, rivaroxaban concentrations were associated with clinical outcome, and the interindividual coefficient of variability was approximately 30% to 40%.¹⁸ This variability stems from the interplay between non-genetic factors, particularly renal insufficiency, advanced age, low body weight, and hepatic impairment, with genetic variations affecting drug transporters and metabolic enzymes.

Pharmacogenetics elucidates the genetic determinants of interindividual variability in drug responses, providing a scientific framework for personalized dosing optimization and adverse event mitigation in antithrombotic therapies. Despite robust mechanistic associations between genetic polymorphisms and rivaroxaban PK, clinical translation remains hindered by heterogeneous evidence and insufficient outcome-driven guidelines.

This review synthesizes evidence from 12 studies based on the following inclusion criteria to evaluate genotype-dependent pharmacokinetic variability and bleeding risk stratification: 1) patients undergoing treatment with rivaroxaban who are of any race, sex or age; 2) reported at least one of the different SNPs on kinetic parameters for subjects or a certain clinical outcome. The exclusion criteria were duplicate publications, literature published in languages other than English, abstracts without essential details and unqualified data. The main characteristics of the included studies are summarized in Tables 1 and 2.

Table 1 Main Characteristics of the Included Studies

Table 2 Visual Representation of Investigated SNP and Outcomes in Each Study

Pharmacokinetic and Pharmacodynamic Profile

Rivaroxaban exhibits high oral bioavailability (80–100%) with dose-dependent food effects; administration with food enhances the bioavailability of the 15/20 mg doses to >80%, while the 10 mg dose remains unaffected. Peak plasma concentrations (C_max) were achieved within 2–4 h of dosing. The drug was highly protein-bound (92–95%, primarily albumin), resulting in a moderate volume of distribution (~50 L).

The elimination of rivaroxaban proceeds via a dual pathway, including the metabolic degradation of the drug and renal elimination of the unchanged drug, see Figure 1. Approximately two-thirds of ingested rivaroxaban is subjected to metabolic degradation by cytochrome P450 enzymes (CYP3A4/5:50% of the total metabolism; CYP2J2:14%) and non-CYP-mediated mechanisms. No pharmacologically active metabolites were detected in the plasma. Renal elimination accounts for one-third of the dose, comprising glomerular filtration (passive) and active tubular secretion mediated by efflux transporter P-glycoprotein (P-gp, encoded by ABCB1) and breast cancer resistance protein (BCRP, encoded by ABCG2).

Figure 1 Proteins and genes associated with the transport and metabolism of rivaroxaban. Abbreviations: P-gp, P-glycoprotein; BCRP, breast cancer resistance protein; CYP450, Cytochrome P450. Notes: ABCB1, CYP3A4 and other italics in the figure and text refer to genes encoding related proteins.ABCB1 loci include rs1045642 (n=12 studies), rs4148738 (n=9), rs1128503 (n=8) and rs2032582 (n=7),rs4728709, rs3789243, rs3213619 (n=1); ABCG2 loci include rs2231142 (n=4), rs2231137 (n=2), rs3114018 (n=1), rs2622604 (n=1), rs1481012 (n=1); ABCA6 locus include rs7212506 (n=1); CYP3A4 loci include rs35599367 (n=2), rs2242480 (n=2), rs4646437 (n=2), rs12333983 (n=1); CYP3A5 locus include rs776746 (n=7), rs15524 (n=1), rs4646450 (n=1); CYP2J2 locus include rs890293 (n=3); CYP2C19 loci include rs4244285 (n=2),rs12248560 (n=2); AKR7A3 loci include rs1738023 (n=1),rs1738025 (n=1).

Pharmacogenetic

The Effects of ABCB1 Polymorphisms

The humans ABCB1 gene³¹ is adjacently located on chromosome 7q21, encoding a drug transporter (P-gp).³² A total of 12 studies investigated ABCB1 genetic variants, with predominant focus on four clinically relevant loci: rs1045642 (n=12 studies), rs4148738 (n=9), rs1128503 (n=8) and rs2032582 (n=7). Secondary loci (rs4728709, rs3789243 and rs3213619) were examined in a single study. There was substantial methodological heterogeneity across the studies, including statistical approaches, bleeding outcome definitions, and pharmacokinetic measurement protocols.

Regarding the types of site mutations, the reported studies did not use a uniform notation system. Therefore, in subsequent descriptions, the same site mutation may be represented in two different ways (Table S1). For example, the site rs1045642 is denoted as c.3435T>C when using the cDNA sequence as a reference and as g.87138645A>G when using the genomic sequence as a reference. However, this does not affect the interpretation of the study results.

rs1045642 (g.87138645A>G)

Among the 12 studies investigating rs1045642 in this review, two yielded discordant ethnicity-specific findings. In Asian populations, heterozygous (AG) and homozygous (GG) mutant genotypes correlated with reductions in rivaroxaban peak plasma concentrations (C_max) compared to wild-type (AA) carriers (P<0.05),¹⁹ whereas no such association was observed in Caucasian cohorts (heterozygous for mutation, β=−8.53, CI: (−32.11, 15.04); homozygous for mutation, β=−12.17, CI: (−32.48 to 8.14).²⁰ This discrepancy likely reflects ethnic disparities in both allele frequencies and baseline drug exposure profiles. Data from the NCBI database’s global population study shows that the gene frequency distribution in European and American populations is A=0.521607, G=0.478393, while in East Asian populations, the gene frequency distribution is A=0.3838, G=0.6162. Meanwhile, research shows that Japanese populations exhibit a 20–40% higher rivaroxaban area under the curve (AUC) than Caucasians.³³ Despite these C_max variations, five studies consistently reported null associations between rs1045642 and trough concentration (C_min) or dose-adjusted residuals across diverse populations (P>0.10).^21–25

Bleeding risk correlations remain contentious owing to the methodological heterogeneity. Sennesael et al²⁶ first identified a non-significant trend toward bleeding in variant allele carriers, although this was limited by a small cohort (n=10). Subsequent studies have reported contradictory findings. Sychev et al²³ reported an increased risk of clinically relevant non-major bleeding (CRNMB) in TT versus CC genotypes (29.3% vs 4.5%, P= 0.021), while Kim et al²⁷ demonstrated a confounder-adjusted elevated major bleeding risk in AA genotype carriers (Model I: 3.243 (1.371–7.671); Model II: 3.167 (1.349–7.436). Conversely, Wang et al²² observed no association among Mongolian descendants (P=0.9107), which potentially reflects ancestry-specific penetrance. Other studies have also failed to observe any correlation.^28,29

Critical limitations include inconsistent outcome definitions, variable adjustments for renal function and P-gp inhibitors, and population stratification effects. These unresolved discrepancies underscore the necessity for standardized phenotyping protocols and ancestry-stratified analyses in future pharmacogenomic investigations.

rs4148738 (g.87163049C>T)

The ABCB1 rs4148738 polymorphism exhibited conflicting pharmacogenetic associations across nine observational studies. Five studies found no significant correlation between the rs4148738 genotypes and C_min in either Asian or Caucasian cohorts.^{19,22,23,26,30} However, subsequent studies revealed population-specific effects: Sychev et al demonstrated that Caucasian patients with the CT genotype had higher residual drug concentrations than CC wild-type carriers (P = 0.039),²⁴ while Wu Tingting et al observed that Asian TT homozygotes showed lower dose-adjusted C_min (C_min/D) compared to CC genotypes [median (IQR) TT: 0.66 (0.31,1.20) vs CC: 1.22 (0.66,3.47); Adjusted P = 0.033]²⁵ It reported that no correlation was observed between peak concentration and this locus.²⁹

Regarding bleeding events, one study²³ reported that the ABCB1 rs4148738 mutant genotype (CT/TT) was associated with an increased incidence of bleeding events compared to the wild-type genotype (CC) (TT vs CC 39.3% vs 8.1%, P = 0.008; TT vs CT 39.3% vs 14.3%, P = 0.002), which is supported by other studies^19,28 In contrast, efficacy assessments focusing on anticoagulation intensity demonstrated no significant alterations in PT levels across the rs4148738 genotypes (P=0.640488).

rs1128503 (g.87179601A>G)

This systematic review incorporated eight studies investigating the ABCB1 rs1128503 polymorphism. Wang et al demonstrated that Mongolian ethnic carriers of the wild-type genotype exhibited a significantly higher C_min than those with mutant genotypes [TT: 33.80 (19.00, 51.90) vs CC: 29.10 (15.61, 57.22), P=0.0421].²² Ain N U et al reported that the heterozygous and homozygous mutant showed lower peak concentrations as compared to the wild-type genotype [AA: 12.06±7.31; AG: 6.51 ±4.15; GG: 5.66±2.42].¹⁹ However, subsequent investigations in multi-ethnic cohorts, including Asian populations^21,25 and European/American populations,²⁰ failed to establish statistically significant associations between this genetic variant and pharmacokinetic parameters.

Notably, while Wang et al identified genotype-dependent concentration variations, no direct genetic correlation with bleeding events was observed (P=0.8062),²² suggesting potential synergistic effects of other pharmacogenetic factors or clinical variables on hemorrhagic outcomes. These findings were corroborated by independent studies by Wu et al,^25–28 which collectively demonstrated the absence of significant genotype-bleeding associations.

Other Loci

Wu et al first identified a clinically significant association between the ABCB1 rs4728709 polymorphism and C_min/D, with heterozygous GA genotype carriers demonstrating significantly reduced C_min/D values compared to wild-type GG individuals (GA vs GG, P = 0.032) [median (IQR): GA: 0.40 (0.30,0.88) vs GG: 1.06 (0.50,1.81); Adjusted P = 0.032].²⁵ This seminal finding provides novel insights into the genetic determinants of rivaroxaban pharmacokinetic variability of rivaroxaban. Notably, complementary evidence from dexamethasone studies revealed that the rs4728709 variant accelerates drug clearance of this P-gp substrate,³⁴ resulting in corresponding reductions in systemic exposure, a mechanistic parallel consistent with Wu’s pharmacometric observations.

Divergent observations regarding ABCB1 rs2032582 exist in the current literature. Zhang et al reported an elevated dose-adjusted C_max (C_max/D) in G allele carriers (P = 0.025, FDR = 0.042)²⁹ and Ain N U et al reported mutant genotypes showed lower peak concentrations as compared to the wild-type genotype,¹⁹ whereas Lenoir C et al found no significant influence of this polymorphism on C_max²⁰ and Nakagawa et al found no significant influence of this polymorphism on C_min/D,²¹ necessitating further validation in larger pharmacogenetic cohorts.²⁶ Meanwhile, no correlation was found in the studies of bleeding events related to this locus.^27,28 Exploratory analyses of less-characterized variants (rs3789243 and rs3213619) have thus far failed to establish significant associations with either rivaroxaban pharmacokinetic parameters or hemorrhagic outcomes (P=0.346; P=0.696).²⁷

In summary, the correlation between ABCB1 polymorphisms and rivaroxaban pharmacokinetics and hemorrhagic events has not been harmonized, and there are considerable differences among different ethnic groups. Most of the available data are from Caucasians, and very little data are available for the Chinese population. Therefore, further research is needed to explore the correlation between ABCB1 gene polymorphisms and P-gp expression, especially in multi-tissue, multi-site, pharmacokinetics and pharmacodynamics, a more systematic study should be carried out on the Chinese population, and further research on the effects of ABCB1 gene polymorphisms and P- gp expression and functional differences on drug disposal, so as to provide a theoretical and practical basis for clinical rational drug use.

The Effects of ABCG2 Polymorphisms

ABCG2 encodes BCRP, an ATP-binding cassette efflux transporter that mediates the cellular extrusion of various substrates including rivaroxaban.⁶ BCRP is predominantly expressed on the apical membrane of intestinal epithelial cells, where it modulates drug bioavailability by actively transporting substrates from the enterocytes into the intestinal lumen. The transporter is also functionally expressed in critical pharmacological barriers and excretory organs, including the blood-brain barrier (limiting central nervous system penetration of substrates), hepatocytes (facilitating hepatobiliary elimination), and renal proximal tubules (mediating active secretion into urine).

Our systematic review included four pharmacogenetic studies that investigated clinically relevant ABCG2 polymorphisms: rs2231142 (n=4), rs2231137 (n=2), rs3114018 (n=1), rs2622604 (n=1), and rs1481012 (n=1).

Regarding the pharmacokinetic outcomes, we have not identified any ABCG2 gene polymorphisms associated with the pharmacokinetics of rivaroxaban. Although ABCG2 rs2231137 and rs2231142 are missense mutations, four studies consistently demonstrated that these loci are not associated with the pharmacokinetic parameters (C_trough/D) of rivaroxaban.^21,25,27,28 Regarding bleeding events, the finding that A carriers of rs3114018 were associated with bleeding complications was supported by Kim (A allele carriers 26.8% vs CC genotype carriers 15.6%, P=0.020).²⁷ Although the SNP is located in an intron non-expressed region, it is possible that those located in the control region of the ABCG2 gene, both in intron and promoter sequences, could affect RNA splicing and thus interfere with the expression/function of ABCG2 proteins, or could result in modified substrate selectivity.³⁵ Kim et al found no significant ABCG2 variants (rs2622604/rs1481012) beyond rs3114018.²⁷ As ABCB1 (P-gp) and ABCG2 (BCRP) exhibit overlapping substrate specificity and synergistic effects, further clinical studies are required to clarify their combined effects on rivaroxaban pharmacokinetic, efficacy and bleeding risk.³²

The Effects of CYP3A4/5 Polymorphisms

Genetic polymorphisms in cytochrome P450 (CYP) enzymes, particularly the CYP3A4/5 isoforms, have been implicated in interindividual variability in drug metabolism.³⁶ Similarly, studies have also demonstrated a strong correlation between CYP3A family activity and rivaroxaban metabolism.³⁷ Therefore, it stands to reason that genetic polymorphisms in the CYP3A family would be associated with the clinical outcomes of rivaroxaban. However, the results of many recent studies contradict this assumption. This review analyzed seven studies that investigated CYP3A4/5 variants (rs35599367, rs2242480, rs4646437, rs12333983 in CYP3A4 and 3* rs776746, rs15524, rs4646450 in CYP3A5).

No significant associations were observed between CYP3A5 3* (rs776746) and rivaroxaban exposure^21,24,30 or between CYP3A4 (rs35599367) and C_min in atrial fibrillation patients, which aligns with previous research findings.^21,24,25,30 Wu et al also did not find that CYP3A4 gene polymorphisms (rs2242480 and rs4646437) had a significant effect on the C_min/D of rivaroxaban. No correlation was observed between CYP3A4/5 (including the above sites and rs12333983, rs15524, rs4646450) and hemorrhage events.^25,27–29 The current data do not support definitive correlations between CYP3A4/5 genetic variants and the PK or efficacy of rivaroxaban.

The Effects of CYP2J2 Polymorphisms

CYP2J2 accounts for approximately 14% of the total clearance of rivaroxaban, which is comparable to the contribution of CYP3A.³⁸ Notably, the catalytic efficiency of CYP2J2 is higher than that of CYP3A4 in vitro. The intrinsic clearance of rivaroxaban catalyzed by CYP2J2 is nearly 39 times greater than that catalyzed by CYP3A4.³⁹ Among the three studies that analyzed CYP2J2 genetic variants in this review, all focused exclusively on the *7 allele (rs890293).

It was noted that there was no correlation between CYP2J2*7 (rs890293) and C_min/D (P=0.331) or C_max/D (P=0.445) in two studies,^21,29 in terms of safety, Zhang et al and Campos et al also did not observe a significant association between CYP2J2*7 (rs890293) and bleeding events (P=0.999).^28,29 CYP2J2 activity is affected by various polymorphisms (such as CYP2J2*2, *3, *4, *6, *8, and *10).³⁹ While current evidence remains inconclusive regarding CYP2J2-mediated pharmacodynamic interactions with rivaroxaban, the polymorphic landscape of the gene warrants further investigation to elucidate the potential subpopulation-specific effects obscured by phenotypic heterogeneity.

The Effects of CYP2C19 Polymorphisms

Although CYP2C19 contributes to rivaroxaban metabolism, pharmacogenomic analyses of two key polymorphisms (CYP2C19*2 [rs4244285] and *17 [rs12248560]) revealed no significant pharmacokinetic associations. A study of patients with atrial fibrillation and acute coronary syndrome demonstrated that neither CYP2C19 *2 nor *17 polymorphism significantly influenced the C_min of rivaroxaban.³⁰ Complementary findings from a separate investigation showed that the CYP2C19 *17 variant was not associated with multiple pharmacokinetic parameters (C_min, C_min/D, C_max/D) and bleeding outcomes.²⁹

Notably, current evidence does not support the clinically relevant CYP450 genetic determinants of rivaroxaban pharmacokinetics. This apparent paradox, in which CYP450 enzymes mediate metabolic clearance yet lacks identified genetic modifiers, may reflect methodological limitations in existing studies. Key constraints include underpowered sample sizes, heterogeneous patient populations, and insufficient characterization of rare CYP450 variants. Systematic pharmacogenomic investigations employing standardized pharmacokinetic phenotyping and multi-ethnic cohorts are required to resolve potential population-specific effects and to elucidate the complex interplay between genetic variation and rivaroxaban metabolism.

Other Genes

The metabolism of rivaroxaban is mainly dependent on the CYP subtypes described above as well as on non-CYP components. AKR7A3 is a member of the AKR family and is involved in exogenous drug metabolism. Considering the complexity of rivaroxaban metabolism, there is a study based on whole-exome sequencing to explore candidate genes associated with the potential bleeding risk of rivaroxaban. Zhao et al found that compared with heterozygous variants and unmutated genotypes, homozygous ABCA6 rs7212506 and AKR7A3 rs1738023/rs1738025 were susceptible sites for bleeding events with rivaroxaban.⁴⁰ ABCA6 is a member of the ATP-binding cassette transporter superfamily that transports both extracellular and intracellular substrates, including drugs and metabolites. They predicted that the AKR7A3 homozygous variant could block the normal metabolism of rivaroxaban, causing systemic accumulation of the active parent drug, while speculating that the ABCA6 homozygous variant might disrupt transmembrane stability and perturb the ABC transporter signaling pathway, thereby altering the disposition of rivaroxaban.

A study⁴¹ exploring the genetic background of DOACs-associated hemorrhagic events found that specific haplotypes (eg, AAAGAGCT and AGAG) were significantly more frequent in hemorrhagic patients than in controls (P<0.05), although no significant differences in any of the examined single-nucleotide variants (SNVs) were detected in hemorrhagic patients. This suggests that haplotype analysis can reveal genetic markers significantly associated with DOACs-related adverse events (AEs) that may not be detected in individual SNV analyses.

Summary

This systematic review synthesizes current evidence on the pharmacokinetic profile of rivaroxaban (encompassing absorption, distribution, metabolism, and excretion characteristics) and pharmacogenetic determinants. Through a critical appraisal of methodological approaches and clinical implications across the included studies, we characterized population-specific variations in drug disposition and evaluated putative genetic biomarkers. Our analysis established an evidence-based framework to (1) identify knowledge gaps in CYP450-mediated metabolic pathways and transporter interactions, (2) quantify the clinical validity of reported genotype-phenotype associations, and (3) assess the evidentiary threshold for implementing genetic testing to optimize rivaroxaban dosing strategies.

Known factors that increase the risk of bleeding during rivaroxaban treatment include renal and hepatic impairments, concomitant therapy with interacting drugs, low body weight, and advanced age.⁴² While pharmacogenetic profiling holds the theoretical potential to guide personalized anticoagulant selection by minimizing interindividual variability in drug exposure, thereby balancing thromboembolic prevention against bleeding risk, current evidence remains insufficient to establish clinically actionable genotype-phenotype correlations. The 25 genetic loci included in the 12 genetic polymorphism studies in this paper include well-known, extensively studied loci as well as newly discovered loci reported in individual studies that are associated with individual differences in rivaroxaban. However, as mentioned earlier, when studying ABCB1 rs1045642, previous reports have shown that the correlation with Cmax varies across different ethnic groups. While many studies have not observed any association with Cmin. First, the data are derived from studies with different statistical methods, results, definitions, and measurements, which may lead to inconsistencies in the results. However, inconsistencies do not mean that such studies are meaningless; they merely highlight the complexity and extremity of clinical data, further emphasising the importance of real-world validation.

In summary, pharmacogenomic monitoring and bleeding risk assessment prior to rivaroxaban administration may help optimise its efficacy and safety in patients. However, whether clinically actionable genotype-phenotype correlations can be established requires larger-scale, robust, global, multicentre clinical trials to validate the potential genetic loci identified by the test and a large-scale data repository to provide the foundation for personalised treatment.

Funding

This study was supported by the Jinhua Central Hospital Youth Research Foundation (Grant No: JY2022-2-03) and Key Science and Technology Plan Projects of Jinhua City (Grant No: 2023-3-122).

Disclosure

The authors report no conflicts of interest in this work.

References

1. Perzborn E, Strassburger J, Wilmen A, et al. In vitro and in vivo studies of the novel antithrombotic agent BAY 59-7939–an oral, direct factor xa inhibitor. J Thromb Haemost. 2005;3:514–521. doi:10.1111/j.1538-7836.2005.01166.x

2. Laux V, Perzborn E, Kubitza D, Misselwitz F. Preclinical and clinical characteristics of rivaroxaban: a novel, oral, direct factor Xa inhibitor. Semin Thromb Hemost. 2007;33:515–523. doi:10.1055/s-2007-982083

3. Gerotziafas GT, Elalamy I, Depasse F, Perzborn E, Samama MM. In vitro inhibition of thrombin generation, after tissue factor pathway activation, by the oral, direct factor Xa inhibitor rivaroxaban. J Thromb Haemost. 2007;5:886–888. doi:10.1111/j.1538-7836.2007.02429.x

4. Anderson DR, Dunbar M, Murnaghan J, et al. Aspirin or Rivaroxaban for VTE prophylaxis after hip or knee arthroplasty. N Engl J Med. 2018;378:699–707. doi:10.1056/NEJMoa1712746

5. Young AM, Marshall A, Thirlwall J, et al. Comparison of an oral factor xa inhibitor with low molecular weight heparin in patients with cancer with venous thromboembolism: results of a randomized trial (SELECT-D). J Clin Oncol. 2018;36:2017–2023. doi:10.1200/JCO.2018.78.8034

6. Bansilal S, Bloomgarden Z, Halperin JL, et al. Efficacy and safety of rivaroxaban in patients with diabetes and nonvalvular atrial fibrillation: The rivaroxaban once-daily, oral, direct factor XA inhibition compared with vitamin K antagonism for prevention of stroke and embolism trial in atrial fibrillation (ROCKET AF Trial). Am Heart J. 2015;170:675–682.e8. doi:10.1016/j.ahj.2015.07.006

7. Male C, Lensing AWA, Palumbo JS, et al. Rivaroxaban compared with standard anticoagulants for the treatment of acute venous thromboembolism in children: a randomised, controlled, Phase 3 trial. Lancet Haematol. 2020;7:e18–e27. doi:10.1016/S2352-3026(19)30219-4

8. Young G, Lensing AWA, Monagle P, et al. Rivaroxaban for treatment of pediatric venous thromboembolism. An Einstein-Jr phase 3 dose-exposure-response evaluation. J Thromb Haemost. 2020;18:1672–1685. doi:10.1111/jth.14813

9. Branch KR, Probstfield JL, Eikelboom JW, et al. Rivaroxaban with or without aspirin in patients with heart failure and chronic coronary or peripheral artery disease. Circulation. 2019;140:529–537. doi:10.1161/CIRCULATIONAHA.119.039609

10. Connor P, Sánchez van Kammen M, Lensing AWA, et al. Safety and efficacy of rivaroxaban in pediatric cerebral venous thrombosis (EINSTEIN-Jr CVT). Blood Adv. 2020;4:6250–6258. doi:10.1182/bloodadvances.2020003244

11. Stevens SM, Woller SC, Kreuziger LB, et al. Antithrombotic Therapy for vte disease: second update of the chest guideline and expert panel report. Chest. 2021;160:e545–e608. doi:10.1016/j.chest.2021.07.055

12. Hindricks G, Potpara T, Dagres N, et al. 2020 ESC Guidelines for the diagnosis and management of atrial fibrillation developed in collaboration with the European Association for Cardio-Thoracic Surgery (EACTS): the Task Force for the diagnosis and management of atrial fibrillation of the European Society of Cardiology (ESC) Developed with the special contribution of the European Heart Rhythm Association (EHRA) of the ESC. Eur Heart J. 2021;42:373–498. doi:10.1093/eurheartj/ehaa612

13. Lawton JS, Tamis-Holland JE, Bangalore S, et al. 2021 ACC/AHA/SCAI guideline for coronary artery revascularization: executive summary: a report of the American College of Cardiology/American Heart Association Joint Committee on Clinical Practice Guidelines. Circulation. 2022;145:e4–e17. doi:10.1161/CIR.0000000000001039

14. Panel P, Abramson BL, Al-Omran M, et al. Canadian cardiovascular society 2022 guidelines for peripheral arterial disease. Can J Cardiol. 2022;38:560–587. doi:10.1016/j.cjca.2022.02.029

15. The Top 300 of 2019. Available from: https://clincalc.com/DrugStats/Top300Drugs.aspx. Accessed August 16, 2025.

16. Yu Z, Yu L, Shan C. Trends of ambulatory oral anticoagulant prescription in five major cities of China, 2012-2017. BMC Health Serv Res. 2020;20:209. doi:10.1186/s12913-020-5072-3

17. Miklič M, Mavri A, Vene N, et al. Intra- and inter- individual rivaroxaban concentrations and potential bleeding risk in patients with atrial fibrillation. Eur J Clin Pharmacol. 2019;75:1069–1075. doi:10.1007/s00228-019-02693-2

18. Testa S, Tripodi A, Legnani C, et al. Plasma levels of direct oral anticoagulants in real life patients with atrial fibrillation: results observed in four anticoagulation clinics. Thromb Res. 2016;137:178–183. doi:10.1016/j.thromres.2015.12.001

19. Ain NU, Ali N, Ullah A, Ullah S, Ahmad S. Bleeding events associated with rivaroxaban therapy in naive patients with nonvalvular atrial fibrillation: A longitudinal study from a genetic perspective with INR follow-up. Medicina. 2024;60:1712. doi:10.3390/medicina60101712

20. Lenoir C, Terrier J, Gloor Y, et al. Impact of the genotype and phenotype of CYP3A and P-gp on the apixaban and rivaroxaban exposure in a real-world setting. J Pers Med. 2022;12:526. doi:10.3390/jpm12040526

21. Nakagawa J, Kinjo T, Iizuka M, et al. Impact of gene polymorphisms in drug‐metabolizing enzymes and transporters on trough concentrations of rivaroxaban in patients with atrial fibrillation. Basic Clin Pharma Tox. 2021;128:297–304. doi:10.1111/bcpt.13488

22. Wang Y, Chen M, Chen H, Wang F. Influence of ABCB1 gene polymorphism on rivaroxaban blood concentration and hemorrhagic events in patients with atrial fibrillation. Front Pharmacol. 2021;12:639854. doi:10.3389/fphar.2021.639854

23. Sychev D, Ostroumova O, Cherniaeva M, et al. The influence of ABCB1 (rs1045642 and rs4148738) gene polymorphisms on rivaroxaban pharmacokinetics in patients aged 80 years and older with nonvalvular atrial fibrillation. High Blood Press Cardiovasc Prev. 2022;29:469–480. doi:10.1007/s40292-022-00536-3

24. Sychev DA, Sokolov AV, Reshetko OV, et al. Influence of ABCB1, CYP3A5 and CYP3A4 gene polymorphisms on prothrombin time and the residual equilibrium concentration of rivaroxaban in patients with non-valvular atrial fibrillation in real clinical practice. Pharmacogenet Genomics. 2022;32:301–307. doi:10.1097/FPC.0000000000000483

25. Wu T, Wu S, Li L, et al. The impact of ABCB1, CYP3A4/5 and ABCG2 gene polymorphisms on rivaroxaban trough concentrations and bleeding events in patients with non-valvular atrial fibrillation. Hum Genomics. 2023;17:59. doi:10.1186/s40246-023-00506-3

26. Sennesael AL, Larock A-S, Douxfils J, et al. Rivaroxaban plasma levels in patients admitted for bleeding events: insights from a prospective study. Thromb J. 2018;16:28. doi:10.1186/s12959-018-0183-3

27. Kim H, Song T-J, Yee J, et al. ABCG2 gene polymorphisms may affect the bleeding risk in patients on apixaban and rivaroxaban. Drug Des Devel Ther. 2023;17:2513–2522. doi:10.2147/DDDT.S417096

28. Campos-Staffico AM, Dorsch MP, Barnes GD, et al. Eight pharmacokinetic genetic variants are not associated with the risk of bleeding from direct oral anticoagulants in non-valvular atrial fibrillation patients. Front Pharmacol. 2022;13:1007113. doi:10.3389/fphar.2022.1007113

29. Zhang D, Qin W, Du W, et al. Effect of ABCB1 gene variants on rivaroxaban pharmacokinetic and hemorrhage events occurring in patients with non-valvular atrial fibrillation. Biopharm Drug Dispos. 2022;43:163–171. doi:10.1002/bdd.2327

30. Sychev DA, Baturina OA, Mirzaev KB, et al. CYP2C19*17 may increase the risk of death among patients with an acute coronary syndrome and non-valvular atrial fibrillation who receive clopidogrel and rivaroxaban. Pharmgenomics Pers Med. 2020;13:29–37. doi:10.2147/PGPM.S234910

31. Ishikawa T, Hirano H, Onishi Y, Sakurai A, Tarui S. Functional Evaluation of ABCB1 (P-Glycoprotein) polymorphisms: high-speed screening and structure-activity relationship analyses. Drug Metab. Pharmacokinet. 2004;19:1–14. doi:10.2133/dmpk.19.1

32. Gong IY, Mansell SE, Kim RB. Absence of both MDR1 (ABCB1) and breast cancer resistance protein (ABCG2) transporters significantly alters rivaroxaban disposition and central nervous system entry. Basic Clin Pharmacol Toxicol. 2013;112:164–170. doi:10.1111/bcpt.12005

33. Gong IY, Kim RB. Importance of pharmacokinetic profile and variability as determinants of dose and response to dabigatran, rivaroxaban, and apixaban. Cana J Cardiol. 2013;29:S24–S33. doi:10.1016/j.cjca.2013.04.002

34. Yang JJ, Cheng C, Devidas M, et al. Genome-wide association study identifies germline polymorphisms associated with relapse of childhood acute lymphoblastic leukemia. Blood. 2012;120:4197–4204. doi:10.1182/blood-2012-07-440107

35. Custodio A, Moreno-Rubio J, Aparicio J, et al. Pharmacogenetic predictors of severe peripheral neuropathy in colon cancer patients treated with oxaliplatin-based adjuvant chemotherapy: a GEMCAD group study. Ann Oncol. 2014;25:398–403. doi:10.1093/annonc/mdt546

36. Zanger UM, Schwab M. Cytochrome P450 enzymes in drug metabolism: regulation of gene expression, enzyme activities, and impact of genetic variation. Pharmacol Ther. 2013;138:103–141.

37. Sychev DA, Vardanyan A, Rozhkov A, et al. CYP3A activity and rivaroxaban serum concentrations in russian patients with deep vein thrombosis. Genet Test Mol Bioma. 2018;22:51–54. doi:10.1089/gtmb.2017.0152

38. Foerster KI, Hermann S, Mikus G, Haefeli WE. Drug–drug interactions with direct oral anticoagulants. Clin Pharmacokinet. 2020;59:967–980. doi:10.1007/s40262-020-00879-x

39. Zhao T, Chen Y, Wang D, et al. Identifying the dominant contribution of human cytochrome p450 2j2 to the metabolism of rivaroxaban, an oral anticoagulant. Cardiovasc Drugs Ther. 2022;36:121–129. doi:10.1007/s10557-020-07129-z

40. Zhao M, Zhang Q, Wang X, et al. Non-synonymous alterations in AKR7A3 and ABCA6 correlate with bleeding in aged patients treated with rivaroxaban. Clin Transl Sci. 2022;15:923–929. doi:10.1111/cts.13205

41. Samoš M, Škereňová M, Nosáľ V, et al. ABCB1 gene single nucleotide variants and haplotypes in atrial fibrillation patients experiencing adverse events on direct oral anticoagulation: a whole gene exome sequencing study. J Cardiovasc Pharmacol. 2025;86:50–59. doi:10.1097/FJC.0000000000001695

42. Cross B, Turner RM, Zhang JE, Pirmohamed M. Being precise with anticoagulation to reduce adverse drug reactions: are we there yet? Pharmacogenomics J. 2024;24(7). doi:10.1038/s41397-024-00329-y

Introduction

Colorectal cancer (CRC) is a major cause of cancer deaths globally, being the third most diagnosed and second deadliest cancer.¹ Its poor prognosis is mainly due to treatment resistance. Despite advances in surgery, chemotherapy, and targeted therapies, survival rates for advanced CRC remain low.² Immune checkpoint inhibitors (ICIs) like pembrolizumab and nivolumab show promise for metastatic CRC with mismatch-repair deficiency (dMMR) and high microsatellite instability (MSI-H), offering potential long-term remission for these patients.³ However, treating CRC with low immunogenicity, such as MSI-L/pMMR tumors, remains challenging. Current research aims to refine patient selection and develop new combination strategies to boost the effectiveness of immune checkpoint inhibitors.⁴ Consequently, more studies are urgently needed to identify biomarkers that predict the success of targeted and immune therapies in CRC and to explore their underlying mechanisms.

Neddylation is a post-translational modification where NEDD8 is attached to specific proteins via a three-step enzymatic process involving E1, E2, and E3 enzymes.⁵ It targets cullin and non-cullin proteins, activating cullin RING ligases and affecting the stability and function of non-cullin proteins.⁶ This pathway is often overactive in cancers.^7–9 Inhibiting neddylation with the NAE inhibitor MLN4924 (pevonedistat) shows strong anticancer effects by promoting cell death^10–13 and enhancing cancer cell sensitivity to various treatments.^14–18

Beyond its impact on cancer cells, neddylation significantly influences the tumor microenvironment (TME) by modulating the activities of immune cells, such as macrophages, dendritic cells (DCs), and T cells, as well as cancer-associated fibroblasts (CAFs), cancer-associated endothelial cells (CAEs), and various other factors.^19,20 The inhibition of neddylation results in the upregulation of PD-L1 expression and the suppression of cancer-associated immunity,²¹ highlighting neddylation as a promising target for cancer therapy. Moreover, neddylation-related patterns and scores offer potential for distinguishing TME characteristics, predicting patient prognosis, and guiding personalized treatment strategies in oncology.^22,23 However, the specific effects of neddylation and neddylation-related genes (NRGs) on CRC remain not well-studied.

Bulk RNA sequencing (RNA-seq) is extensively utilized in translational research to ascertain the average transcript expression within heterogeneous cell populations.²⁴ In contrast, single-cell RNA sequencing (scRNA-seq), a more recent advancement, evaluates gene expression at the individual cell level, thereby elucidating the distribution of expression across various cell subtypes.²⁵ The application of machine learning techniques enhances the analysis of data derived from both methodologies, thereby refining insights into diseases such as cancer and facilitating their clinical application.^26–28

In this study, we conducted a comprehensive mapping of neddylation-related genes (NRGs) across 16 distinct cell types in colorectal cancer (CRC) at the single-cell resolution. We performed an analysis of functional variations in neddylation patterns and employed ten machine learning algorithms to develop a gene signature, referred to as NRGS. This signature is designed to predict patient prognosis, characterize the tumor microenvironment, and assess responses to immunotherapy and chemotherapy, utilizing bulk RNA sequencing data.

Methods

Data Acquisition

Single-cell transcriptome datasets GSE166555 and GSE139555 were accessed from the Tumor Immune Single-cell Hub (TISCH) database (http://tisch.comp-genomics.org/home/). Bulk RNA-seq data were obtained from the UCSC Xena browser (https://xenabrowser.net/) and the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The TCGA-COADREAD datasets, which include RNA-seq expression matrices, clinical data, and masked annotated somatic mutations, were acquired from the UCSC repository. Additionally, the GSE39582 dataset, comprising 566 colorectal cancer (CRC) samples, and the GSE17538 dataset, containing 238 CRC samples with corresponding clinical information, were downloaded from the GEO database. Neddylation-related genes (NRGs) and gene sets from five pathways (Biocarta, Hallmark, KEGG, Reactome, and Wiki Pathways) were compiled from the Molecular Signatures Database (MsigDB, https://www.gsea-msigdb.org/gsea/msigdb/).

Single-Cell RNA-Seq Analysis

Two single-cell transcriptomic datasets, GSE166555 and GSE139555, were integrated for subsequent analysis utilizing the “Seurat” package,²⁹ which is specifically designed for single-cell RNA sequencing analysis. To address batch effects, regularized negative binomial regression was employed. Non-linear dimensionality reduction was performed using Uniform Manifold Approximation and Projection (UMAP). Cluster-specific biomarkers were identified using the “COSG” package, while differential gene expression for each cell type was determined via the “FindAllMarkers” function. Furthermore, hallmark gene enrichment scores were evaluated using single-sample Gene Set Enrichment Analysis (ssGSEA) from the “GSVA” package, employing hallmark gene sets sourced from the Molecular Signatures Database (MsigDB).

Neddylation Patterns Identified via Unsupervised Clustering

We performed a univariate Cox regression analysis to identify NRGs that significantly impacted overall survival (OS) in CRC (p<0.01), utilizing the “survival”, “survminer”, and “ggplot2” packages. Subsequently, we employed the “ConsensusClusterPlus” package³⁰ to perform unsupervised consensus clustering on the expression profiles of prognostic NRGs. The optimal number of clusters was determined based on the cumulative distribution curve and K-means clustering. The robustness of the clustering results was further validated through principal component analysis (PCA).

Signature of Neddylation-Related Genes (NRGS) Obtained Through Integrative Machine Learning methods

Based on the Univariate Cox regression analysis with a p-value threshold of <0.05, we identified 52 NRGs with significant prognostic value. To develop a robust prognostic NRG signature (NRGS), these biomarkers were subjected to an integrative machine learning analysis. This procedure included various algorithms such as random survival forests (RSF), least absolute shrinkage and selection operator (LASSO), gradient boosting machine (GBM), survival support vector machine (Survival-SVM), supervised principal components (SuperPC), ridge regression, partial least squares regression for Cox (plsRcox), CoxBoost, Stepwise Cox, and elastic network (Enet). Notably, RSF, LASSO, CoxBoost, and Stepwise Cox have capabilities for dimensionality reduction and variable selection, and they were combined with other algorithms to form multiple machine learning algorithm combinations. Similar methodologies have been employed in previous studies.³¹ The TCGA-COADREAD dataset was utilized as the training dataset, while the GSE17538 and GSE39582 datasets served as external validation datasets. The concordance index (C-index) for each model across all validation datasets was computed, and the model exhibiting the highest average C-index was identified as the optimal model. Following the identification of the optimal model, the median risk score derived from the training dataset was employed as the threshold to stratify patients in both the training and validation datasets into high-risk and low-risk groups. Kaplan-Meier survival analysis and the Log rank test were conducted on these two groups using the “survival” and “survminer” R packages. Receiver operating characteristic (ROC) curves were generated to assess the predictive accuracy of the model.

Gene Set Variation Analysis (GSVA)

The hallmark gene sets associated with neddylation and related pathways were obtained from the MsigDB. The enrichment scores of these hallmark genes were assessed using single-sample Gene Set Enrichment Analysis (ssGSEA) via the R package “GSVA.”³²

Integrated Omics Analysis

The “maftools” package³³ was utilized to characterize somatic mutations in genes of CRC patients.

Functional Enrichment Analysis

To explore the potential biological functions of differential expressed genes, Gene Ontology (GO) functions, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathway enrichment analyses were performed using GSEA with R package “clusterProfiler”.³⁴

Immune Infiltration Analysis

The single-sample GSEA (ssGSEA) method was employed with the R package “GSVA” to evaluate immune cell scores. Furthermore, the R package “IOBR”³⁵ which integrates eight algorithms (MCPcounter, EPIC, xCell, CIBERSORT, IPS, quanTIseq, ESTIMATE, and TIMER), was used to estimate the abundance of immune cells across different risk groups. Moreover, we analyzed hub gene expression in immune and non-immune cells within CRC single-cell datasets using the TISCH database.

Immunotherapeutic Response Prediction and Drug Sensitivity Analysis

To assess the role of NRGS in predicting the benefits of immunotherapy, we employed the Tumor Immune Dysfunction and Exclusion (TIDE) score, accessible via the TIDE website (http://tide.dfci.harvard.edu). Subsequently, we investigated drug susceptibility across the two NRGS-defined risk groups by calculating the half-maximal inhibitory concentration (IC50) values using the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) with the “oncoPredict” R package.³⁶

Statistical Analysis

Data processing, statistical analyses, and visualization were performed using R software (version 4.2.2). The Pearson correlation coefficient was utilized to evaluate the correlation between continuous variables. Chi-square tests were applied for comparisons of categorical variables, while Wilcoxon rank-sum tests or T-tests were used for continuous variables. The “survminer” package was employed to determine optimal cutoff values. Cox regression and Kaplan-Meier analyses were conducted using the “survival” package. Figures were generated with the “ggplot2” package. A p-value of less than 0.05 was considered statistically significant (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).

Results

The Allocation of Various Cell Subsets in CRC According to the scRNA-Seq Atlas

We utilized the scRNA-seq datasets (GSE166555 and GSE139555) and selected 14 samples from CRC patients for further investigation. Rigorous quality control procedures were implemented, resulting in the acquisition of 43,725 cells from the predetermined samples. Through UMAP analysis and hierarchical clustering, we identified 16 distinct cell subgroups, including B cells, CD4+ T cells, CD8+ T cells, endothelial cells, epithelial cells, fibroblasts, dendritic cells (DC), malignant cells, mast cells, monocytes/macrophages (Mono/Macro), myofibroblasts, natural killer (NK) cells, plasma cells, Th17 cells, proliferating T cells (Tprolif), and regulatory T cells (Treg) (Figure 1A). Subsequently, we examined the expression patterns of marker genes across these 16 cell types and presented the top three genes in a heat map (Figure 1B). Using the FindAllMarkers function, we identified differentially expressed genes for each of the 16 cell types. Figure 1C illustrates the top five up-regulated and down-regulated genes for each cell type. The GSVA R package was employed to compute the gene set variation analysis scores for 50 hallmark pathways across the cell types, revealing significant enrichment in endothelial, epithelial, fibroblasts, myofibroblasts, DC, Mono/Macro, and malignant cells (Figure 1D).

Figure 1 Different cell clustering based on scRNA-seq data of CRC and further analysis. (A) Cluster annotation and 16 cell types identification by means of UMAP. (B) Heat map of marker genes for the 16 cell types. (C) Differential expression genes for the 16 cell types. (D) Heat map of the scores of HALLMARK pathways for the 16 cell types. scRNA-seq, single cell RNA sequencing; UMAP, Unified Flowform Approximation and Projection.

The NRGs Status Across Various Cell Subsets as Determined by the scRNA-Seq Atlas

Utilizing the MSigDB database, a total of 250 neddylation-related genes (NRGs) were identified. Figure 2A illustrates the expression patterns of these 250 NRGs across 16 distinct cell types. Gene Set Variation Analysis (GSVA) was employed to compute the neddylation gene set scores for each cell type, categorizing them into high-score and low-score groups. The results were visualized using Uniform Manifold Approximation and Projection (UMAP) plots (Figure 2B). Notably, cell types such as Mono.Macro, DC, NK, malignant, endothelial, fibroblasts, myofibroblasts, and Tprolif exhibited relatively higher neddylation scores (Figure 2C). Figure 2D presents the composition of cell types within the two groups and across individual samples. Significant differences in cell type composition were observed between the high-score and low-score groups. The high-score group demonstrated increased proportions of malignant, Mono.Macro, fibroblasts, endothelial, myofibroblasts, DC, NK, Tprolif, and Treg cells, whereas the low-score group contained relatively higher proportions of CD4 T, B, and CD8 T cells. These findings were further elucidated through linear plots and UMAP plots, offering a comprehensive visualization of cell type distribution (Figure 2E and F).

Figure 2 The status of NRG in different cell subsets based on scRNA-seq atlas. (A) The expression patterns of 250 NRG in 16 cell types. UMAP plots (B) and Box plots (C) showed differential NRG scores in 16 cell types. (D) Box plots showed the composition of the 16 cell types in high and low NRG score groups and in each sample. Linear plots (E) and UMAP plots (F) providing a comprehensive visualization of the 16 cell types distribution in two NRG score groups.

The HALLMARK Pathways Between High-and Low-Score Groups Based on scRNA-Seq Atlas

We subsequently conducted an analysis of the differential enrichment of hallmark pathways between groups with high and low scores. The findings indicated that nearly all pathways were significantly enriched in the high-score group, including the MYC targets pathway, MTORC1 signaling pathway, TGF-beta signaling pathway, KRAS signaling pathway, PI3K-AKT-MTOR signaling pathway, DNA repair, and oxidative phosphorylation pathways (Figure 3A). Furthermore, we examined the correlation between the neddylation score and hallmark pathways across various cell types. The results demonstrated that the neddylation score was positively correlated with almost all pathways in 16 cell types, particularly in endothelial, epithelial, fibroblast, dendritic, malignant, mast, monocyte/macrophage, myofibroblast, plasma, and proliferating T cells (Figure 3B). The most significant pathways included the MYC targets pathway, MTORC1 signaling pathway, PI3K-AKT-MTOR signaling pathway, DNA repair, and oxidative phosphorylation pathways.

Figure 3 The HALLMARK pathways enrichment analysis in two NRG risk score groups based on scRNA-seq atlas. (A) The HALLMARK pathway enrichment between high and low NRG score groups. (B) Correlation of NRG score with HALLMARK pathways in 16 cell types. * p<0.05, ** p<0.01, *** p<0.001, **** p<0.0001.

Discovery of Various Neddylation Patterns in CRC Using Bulk-RNA Sequencing Data

Subsequently, we conducted a univariate Cox regression analysis on 250 NRGs in CRC samples, which led to the identification of 27 NRGs with significant prognostic value (p<0.01) that exhibited inter-correlation (Figure 4A). Utilizing a consensus clustering algorithm, optimal clustering was achieved at k=2 (Figure 4B), allowing us to categorize the CRC samples into two distinct neddylation patterns, designated as Cluster A (n=701) and Cluster B (n=455). Notably, patients in Cluster B demonstrated poorer overall survival compared to those in Cluster A (Figure 4C). Furthermore, there were significant differences in the expression levels of all prognostic NRGs, with the exception of COPS7A, between the two clusters (Figure 4D). Specifically, Cluster A exhibited higher expression levels of CISH, CCNF, BIRC5, PSME3, PSMB8, FBXW9, PSMA5, BRCA1, PSMA7, PSME1, DCAF7, PSMD3, PSMB10, and PSMB2, whereas Cluster B showed elevated expression of SOCS5, KLHL20, CCDC8, FBXL7, WSB1, FEM1B, COPS8, FBXO30, RNF7, COPS2, FBXO15, and FBXO32. Subsequently, we generated a heat map to visualize the clinicopathological factors and the expression of the 27 prognostic NRGs across the two clusters in CRC patients (Figure 4E). These findings support the hypothesis that distinct neddylation patterns are indicative of varying prognoses and clinical characteristics.

Figure 4 Identification of two neddylation patterns of CRC based on bulk-RNA sequencing data. (A) The correlation network of 27 NRG with significant prognostic value in CRC identified by univariate Cox analysis with p-value <0.01 as the cut off. In the right half, purple denotes a prognostic risk factor, and green signifies a protective factor. Circle size reflects the P value magnitude. A line indicates a correlation between related genes with p < 0.05. (B) The optimal number of clusters based on k=2. (C) Kaplan-Meier analysis of prognosis between the two patterns. (D and E) Differences in clinical characteristics and prognostic NRG expression between the two patterns. * p<0.05, *** p<0.001, ns, not significant.

The Enrichment of Pathways Between the Two Neddylation Patterns

We subsequently acquired five pathway gene sets from the MSigDB database, specifically Biocarta, Hallmark, KEGG, Reactome, and WikiPathways. The R package “GSVA” was employed to score these pathways, and the R package “pheatmap” was utilized to generate a heatmap for comparing the two groups. Within the Biocarta pathway, Cluster A exhibited enrichment in pathways related to MCM, P27, G2, cell cycle, mitochondrial function, P53, MHC, Fas, and TNFR1. Conversely, Cluster B demonstrated enrichment in pathways associated with NK cells, CD40, TCRA, IL-6, ERK5, BCR, IL-10, IL-4, IL-3, MET, TCR, CXCR4, IL1R, IGF1, IL-7, PTEN, ECM, ALK, B lymphocytes, CCR5, TGF-β, LYM, and monocytes. Regarding the Hallmark pathway, Cluster A was enriched in MYC targets, G2M checkpoint, oxidative phosphorylation, DNA repair, and MTORC1 signaling. In contrast, Cluster B was enriched in IL6-JAK-STAT3 signaling, apoptosis, IL2-STAT5 signaling, NOTCH signaling, hypoxia, inflammatory response, TNFA signaling via NFKB, KRAS signaling, TGF-β signaling, angiogenesis, and EMT. In the context of KEGG pathways, Cluster A demonstrated significant enrichment in pathways related to DNA replication, mismatch repair, base excision repair, the cell cycle, and oxidative phosphorylation. Conversely, Cluster B exhibited enrichment in pathways associated with B cell receptor signaling, chemokine signaling, cytokine-cytokine receptor interactions, MAPK signaling, adherens junctions, TGF-β signaling, focal adhesion, and ECM receptor interactions. Regarding Reactome pathways, Cluster A was enriched in processes such as DNA replication, mismatch repair, G2-M checkpoints, base excision repair, cell cycle checkpoints, G2-M DNA damage checkpoints, stabilization of P53, DNA repair, cellular response to hypoxia, regulation of TP53 activity, apoptosis, and MTORC1-mediated signaling. In contrast, Cluster B was enriched in pathways involving IL-3, IL-5, and GM-CSF signaling, IL-27 signaling, IL-21 signaling, PI3K-AKT signaling in cancer, RET signaling, IL-15 signaling, AKT signaling in cancer, and the IL-6 signaling pathway. In the WikiPathway analysis, Cluster A exhibited significant enrichment in pathways related to DNA replication, DNA mismatch repair, base excision repair, metabolic reprogramming in colon cancer, the cell cycle, and ATM signaling. Conversely, Cluster B demonstrated enrichment in pathways associated with TGF-β signaling, B cell receptor signaling, T cell receptor signaling, NOTCH signaling, chemokine signaling, IL-3 signaling, Wnt signaling, IL-4 signaling, MAPK signaling, RAS signaling, IL-10 anti-inflammatory signaling, PI3K-AKT signaling, focal adhesion, TGF-β receptor signaling, CCL18 signaling, and angiogenesis. These findings indicate distinct immune activities and proliferative statuses between the two risk groups, potentially explaining the observed differences in survival rates (Figure 5).

Figure 5 Pathways enrichment analysis in two neddylation patterns. The pathways enrichment of Biocarta (A), Hallmark (B), KEGG (C), Reactome (D) and wiki pathway (E) between the two patterns.

Analyses of Somatic Mutation in Two Neddylation Patterns

We subsequently examined the differences in genomic mutations associated with the two neddylation patterns. The somatic mutation landscape of each CRC sample within the cluster A and B groups was visualized using a waterfall diagram. It was observed that the top ten mutant genes in the two subgroups were largely consistent (Figure 6A and B). In the cluster A group, the genes with the highest mutation frequencies were APC (61%), TP53 (57%), TTN (50%), KRAS (46%), and MUC16 (28%). Similarly, in the cluster B group, the most frequently mutated genes were APC (60%), TP53 (55%), KRAS (42%), TTN (40%), and SYNE1 (27%). Regarding the classification of variations, missense mutations, nonsense mutations, and multi-hit mutations were the three most prevalent types across all mutations. We further analyzed the differential mutant genes between the two groups. Figure 6C illustrates the top ten genes with higher mutation frequencies in the cluster A group compared to the cluster B group, namely RIMS1, SPEN, FAT1, FSIP2, YTHDC1, ZNF516, MXRA5, ASTN2, PCNT, and PCDHGA6. Additionally, tumor mutational burden (TMB) was compared between the two groups, and no significant difference was observed (data not shown).

Figure 6 Somatic mutation analysis in two neddylation patterns. (A) The top 10 mutated genes in the cluster A group. (B) The top 10 mutated genes in the cluster B group. (C) Differences in mutated genes between the two patterns. ** p<0.01, *** p<0.001.

The Immune Landscapes Differ Between the Two Neddylation Patterns

Principal Component Analysis (PCA) demonstrated that the various subtypes were heterogeneous and potentially indicative of differing levels of neddylation modifications (Figure 7A). Furthermore, these findings were corroborated using single-sample Gene Set Enrichment Analysis (ssGSEA), revealing that the cluster A group typically exhibited a lower abundance of immune cells, including activated B cells, activated dendritic cells, macrophages, mast cells, Myeloid-derived suppressor cells (MDSCs), NK cells, NKT cells, regulatory T cells, and type 1 T helper cells (Figure 7B). Subsequently, eight algorithms were employed to investigate the variations in the tumor microenvironment (TME) status across the two neddylation patterns. The analysis indicated that the cluster A group had lower immune, stromal, and ESTIMATE scores, but a higher tumor purity score. This group also showed a relatively low abundance of infiltrated immune cell types (Figure 7C). In summary, both analytical approaches suggested that cluster A exhibited reduced levels of immune infiltration, which may contribute to tumor immune evasion.

Figure 7 Different immune landscapes of two neddylation patterns. (A) PCA in two patterns. (B) Abundance of 23 immune cells based on ssGSEA algorithm in two patterns. (C) Eight immune infiltration algorithms exhibit different numbers of immune cells between the two patterns. * p<0.05, ** p<0.01, *** p<0.001, ns, not significant.

Development and Verification of Prognostic Signatures Associated with Neddylation-Related Genes (NRGS)

Utilizing a Univariate Cox regression analysis with a p-value threshold of <0.05, we identified 52 NRGs with significant prognostic implications, comprising 24 risk-associated genes and 28 protective genes. These 52 potential prognostic biomarkers were subsequently subjected to an integrative machine-learning procedure involving 10 previously mentioned methods. We computed the average concordance index (C-index) for each algorithm across all cohorts and selected the StepCox[both]+RSF algorithm, which exhibited the highest average C-index of 0.71, as the optimal model (Figure 8A). Subsequently, we calculated the risk score for each sample within the cohort based on the expression profiles of the 26 NRGs included in the StepCox[both]+RSF model, namely PLAB2, DCAF7, PSMA5, ASB9, COPS7A, FBXO15, COPS8, FBXL16, RNF7, CCDC8, PSMD12, FBXL19, BIRC5, BTRC, PSME1, TRIM40, FBXL7, DCAF5, CCNF, CUL4B, NAE1, KLHL21, SOCS3, DCUN1D3, FBXO32, and OBSL1 (Figure 8B). Kaplan-Meier survival analysis for overall survival (OS) demonstrated that the high-risk score group exhibited significantly poorer survival outcomes in the training cohort (TCGA-CRC) (Figure 8C). Importantly, the area under the receiver operating characteristic (ROC) curve (AUC) for 1-year, 2-year, and 3-year OS was 0.981, 0.987, and 0.987, respectively, in the training cohort (Figure 8D). Moreover, we showed the changes in AUC at different time points in the training cohort and found that the AUC is all close to 1 (Figure 8E). High-score group also had poorer overall survival in the testing cohort (GSE17538 and GSE39582) and all combined cohorts (Figure 8F–H).

Figure 8 A prognostic NRGS developed by machine learning analysis. (A) The C-index of 107 kinds of prognostic models constructed by 10 machine learning algorithms in training and testing cohort. (B) Error rate curve of random forest tree model. (C) Kaplan-Meier analysis of different NRGS score in training cohort (TCGA-COADREAD). (D) The area under the ROC curve (AUC) for the 1-year, 2-year, and 3-year OS in the training cohort. (E) The changes in AUC value of different time points. Kaplan-Meier analysis of different NRGS score in testing cohort (GSE17538) (F), (GSE39582) (G) and in merge cohort (H).

We also investigated the differences in genomic mutations between the high- and low- risk groups. TP53 (60%), APC (53%), KRAS (52%), TTN (45%) and SYNE1 (29%) were the top 5 genes with the highest mutation frequencies in the high-risk group (Figure S1A), while APC (61%), TP53 (56%), TTN (46%), KRAS (43%) and SYNE1 (27%) in the low risk groups (Figure S1B). Missense mutation, nonsense mutation and multi hit are also the top 3 across all mutation types. Figure S1C showed the top 10 higher mutant genes (ENPEP, DLCK1, DCAF12, GAREML, USP12, SUN2, CD74, LYNX1, ADORA3 and GAGNG3) in high-risk group compared with low-risk group.

We utilized the TISCH database to investigate the distribution of the top ten genes (ASB9, CCDC8, COPS7A, COPS8, DCAF7, FBXL16, FBXO15, PALB2, PSMA5, and RNF7) associated with the RSF model at the single-cell level (Figure S2). Our analysis revealed that the distributions of COPS7A, COPS8, DCAF7, PSMA5, and RNF7 were widespread across most detected cell types, particularly in Tprolif cells.

Correlation Between Clinical Features and NRGS

Furthermore, significant differences were observed in the distribution of high and low-score groups concerning recurrence or metastasis, stage, and tumor type (Figure 9A–D). Subsequently, we categorized all patients based on various clinical characteristics and found that the risk score significantly increased with recurrence or metastasis, advanced stage, and COAD (Figure 9A–D). Notably, as illustrated in Figure 9E and F, both univariate and multivariate Cox regression analyses demonstrated that the NRGS-based risk score served as an independent risk factor for the overall survival rate of CRC.

Figure 9 Correlation between Clinical Features and NRGS. The differences in risk scores across clinical subgroups, including gender (A), recurrence or metastasis (B), stage (C) and tumor type (D). Univariate (E) and multivariate Cox regression analysis (F) illustrated that the NRGS could be used as an independent prognostic factor for CRC patients (P< 0.001).

Functional Enrichment Analysis of NRGS

GSEA was employed to further elucidate the biological functions and signaling pathways associated with NRGS. Initially, we examined the genes correlated with the NRGS-based risk score in CRC. Figure 10A and B illustrate the top 50 genes demonstrating positive and negative correlations with the risk score, respectively. Regarding GO terms, genes related to the NRGS risk score were enriched in processes such as extracellular matrix organization, extracellular structure organization, collagen metabolic process, blood vessel development, angiogenesis, negative regulation of tumor necrosis factor production, cell-substrate adhesion, cell-matrix adhesion, and negative regulation of tumor necrosis factor superfamily cytokine production (Figure 10C). In terms of the KEGG pathways, genes associated with the NRGS risk score were enriched in ECM-receptor interaction, phagosome, focal adhesion, cell adhesion molecules, leukocyte transendothelial migration, protein digestion and absorption, calcium signaling pathway, and proteoglycans in cancer (Figure 10D). For Reactome pathways, NRGS risk score-related genes showed enrichment in extracellular matrix organization, amplification of signal from the kinetochores, mitotic spindle checkpoint, G2/M checkpoints, and DNA strand elongation (Figure 10E).

Figure 10 Functional Enrichment Analysis of NRGS. (A) The top 50 genes positively correlated with NRGS score. (B) The top 50 genes negatively correlated with NRGS score. (C) GO functional enrichment analysis of NRGS score related genes. (D) KEGG pathway enrichment analysis of NRGS score related genes. (E) Reactome pathway enrichment analysis of NRGS score related genes.

The TME Status Differ Between the Two NRGS Risk Groups

We employed eight algorithms to investigate the variations in the TME status between high-risk and low-risk groups. The findings from the ESTIMATE algorithm indicated that stromal, immune, and ESTIMATE scores were significantly elevated in the high-risk group compared to the low-risk group, whereas tumor purity was notably higher in the low-risk group. Additionally, cell populations such as activated dendritic cells (aDC), cancer-associated fibroblasts (CAFs), central memory CD4+ T cells (CD4-Tcm), conventional dendritic cells (cDC), cytotoxic lymphocytes, dendritic cells (DC), endothelial cells, immature dendritic cells (IDC), macrophages, M0 macrophages, M2 macrophages, major histocompatibility complex molecules (MHC), monocytes, mesenchymal stem cells (MSC), myeloid dendritic cells, neutrophils, and CD8+ T cells were more prevalent in the high-risk group. Conversely, cell types such as CD4+ memory T cells, CD4+ naive T cells, immunosuppressive cells (SC), immune checkpoint pathways (IPS), resting natural killer (NK) cells, natural killer T (NKT) cells, plasmacytoid dendritic cells (pDC), T helper 1 (Th1) cells, T helper 2 (Th2) cells, and azurophilic granules (AZ) were more predominant in the low-risk group (Figure 11A).

Figure 11 Correlation between immune microenvironment and NRGS in CRC. (A) The correlation between NRGS and the immune cell infiltration is based on eight algorithms. (B) The expression of chemokines and receptors, interieuknes and receptors, interferons and receptors, and other cytokines in CRC patients with different NRGS scores. * p<0.05, ** p<0.01, *** p<0.001.

Subsequently, we conducted an analysis of the differential expression levels of chemokines and their receptors, interleukins and their receptors, interferons and their receptors, as well as other cytokines between high- and low-risk groups. Regarding chemokines and their receptors, the expression levels of CCL5, CCL8, CCL17, CCL18, CCL21, CCL22, CCR1, CCR2, CCR5, CCR10, CXCL12, CXCL16, and CXCR4 were elevated in the high-risk group, whereas CCL20, CCL28, CCR6, and CXCL11 exhibited higher expression in the low-risk group. In terms of interleukins and their receptors, IL1R1, IL1R2, IL4R, IL6, IL10, IL10RA, IL10RB, IL11, IL17B, IL17D, and IL21R showed increased expression in the high-risk group, while IL1A, IL5, IL12A, and IL26 were more highly expressed in the low-risk group. For interferons and their receptors, IFNAR2 and IFNGR2 were more highly expressed in the high-risk group. Among other cytokines, CSF1, CSF2RB, PDGFC, PDGFD, PDGFRA, PDGFRB, TGFB3, TGFBR1, TGFBR2, TNF, VEGFB, VEGFC, and EPOR exhibited higher expression levels in the high-risk group, with EGF being the only cytokine with elevated expression in the low-risk group (Figure 11B).

Analysis of Predicted Responses to Immunotherapy of NRGS

We assessed the Tumour Immune Dysfunction and Exclusion (TIDE) scores of CRC patients using the TIDE database. Our analysis revealed that the high-risk group exhibited significantly higher TIDE scores (Figure 12A). Bar plots depicting the distribution of responders and non-responders within both high- and low-risk groups indicated that 75% of individuals in the high-risk group were non-responders, while 25% were responders. In contrast, the low-risk group comprised 51% non-responders and 49% responders (Figure 12B). Furthermore, the high-risk group demonstrated an increased propensity for tumour immune exclusion (Figure 12C) and dysfunction (Figure 12D).

Figure 12 Prediction of immunotherapy effects of NRGS in CRC. (A) The TIDE score in high- and low-risk groups. (B) The immunotherapy response rate in high- and low-risk groups. The T-cell exclusion score (C) and T-cell dysfunction score (D) in high- and low-risk groups.

Drug Sensitivity of NRGS

Subsequently, the oncoPredict algorithm was employed to estimate the IC50 values, facilitating the prediction of differential chemotherapy responses between high-risk and low-risk patient groups. A significant variance in the sensitivity to numerous chemotherapeutic agents was observed between these groups. Specifically, the results indicated that patients in the high-risk group exhibited lower sensitivity to oxaliplatin (p<0.0001), 5-Fluorouracil (p<0.01), KRAS G12C inhibitor (p<0.0001), sorafenib (p<0.0001), trametinib (p<0.05), dabrafenib (p<0.01), afatinib (p<0.0001), cisplatin (p<0.05), and crizotinib (p<0.001) compared to those in the low-risk group, as illustrated in Figure 13.

Figure 13 Prediction of drug sensitivity of NRGS in CRC. The Y-axis showed the IC50 value, which was negatively correlated with drug sensitivity. *p<0.05, ** p<0.01, ****p < 0.0001.

Discussion

Patients with advanced metastatic colorectal cancer (mCRC) frequently miss the opportunity for surgical intervention, rendering pharmacological treatment the primary therapeutic strategy. In recent years, immune checkpoint inhibitors (ICIs) have demonstrated remarkable and sustained survival benefits in a limited subset of patients characterized by high immunogenicity, specifically those with microsatellite instability-high (MSI-H) or mismatch repair-deficient (dMMR) tumors.³ Conversely, the majority of colorectal cancer (CRC) patients exhibit low immunogenicity, such as microsatellite instability-low (MSI-L) or proficient mismatch repair (pMMR), and consequently show resistance to ICIs. Furthermore, numerous studies have sought to identify several distinct CRC subtypes to predict the prognosis and efficacy of ICIs treatment, including those classified by the Consensual Molecular Subtype (CMS) classification, as well as mutations in DNA polymerase D1 (POLD1) or DNA polymerase E (POLE), and Lynch syndrome (LS), among others.⁴ Therefore, there remains a critical and unmet need to identify molecular subtypes or biomarkers that can enhance the selection of patients likely to respond to specific treatments. Neddylation, a post-translational modification, is frequently hyperactive in cancer and significantly impacts the tumor microenvironment (TME) by modulating immune cell activities, positioning it as a promising target for cancer therapy.^19,20 However, systematic investigations into neddylation-related genes (NRGs) in CRC are currently lacking. Therefore, we conducted a systematic analysis of NRGs in CRC to elucidate their potential roles in tumor progression, prognosis, the TME and treatment responses. This was achieved through the characterization of the single-cell landscape and the integration of bulk transcriptome data.

Initially, we integrated single-cell data from the GSE166555 and GSE139555 cohorts to analyze cellular heterogeneity in CRC. Our analysis indicated that the group with a high neddylation score exhibited increased proportions of malignant cells, monocytes/macrophages, fibroblasts, endothelial cells, myofibroblasts, DC cells, NK cells, and Tprolif cells, all of which are essential components of the tumor microenvironment (TME). These TME components are instrumental in facilitating in tumor progression, and targeting these cells could be crucial in influencing tumor outcomes. Notably, Tprolif cells exhibited the highest neddylation scores and were characterized by elevated expression levels of immunological exhaustion markers, including PDCD1, HAVCR2, CTLA4, LAG3, and TIGIT.³⁷ Importantly, an increasing body of research underscores the significant role of neddylation in regulating the TME.^19,20 These findings suggest that neddylation may significantly influence tumor-stroma interactions in CRC.³⁸

In a subsequent analysis, we identified 27 neddylation-related genes (NRGs) with significant prognostic value (p<0.01) through univariate Cox regression analysis. Employing a consensus clustering algorithm, we further delineated two distinct neddylation patterns, termed Cluster A and Cluster B. Notably, differential expression was observed in all prognostic NRGs, with the exception of COPS7A, between these two clusters. Cluster B was associated with a poorer prognosis. Single-sample Gene Set Enrichment Analysis (ssGSEA) indicated that Cluster B exhibited a higher abundance of immune cells, including B cells, T cells, and natural killer (NK) cells. However, it also exhibited an increased presence of cells functioning as immunosuppressive regulators in cancer, such as myeloid-derived suppressor cells (MDSCs),³⁹ regulatory T cells (Tregs),⁴⁰ natural killer T (NKT) cells,⁴¹ and tumor-associated macrophages (TAMs).^42–44 Interestingly, Treg cells have been observed to enhance TAM activity, thereby inducing a suppressive TME through the inhibition of CD8+ T cells.⁴⁵ This cellular composition may contribute to a more aggressive cancer phenotype and a less favorable prognosis for patients. Additionally, Cluster B demonstrated a significant enrichment in several signaling pathways, including TGF-β, IL6-JAK-STAT3, PI3K-AKT, NOTCH, MAPK, RAS, and angiogenesis pathways. Extensive research has established that the aberrant activation of these pathways plays a pivotal role in the proliferation, migration, and invasion of CRC, thus altering the TME and contributing to a poorer prognosis.^46–51 Notably, neddylation is crucial in modulating the TME by influencing various signaling pathways, such as TGF-β, PI3K-AKT-mTOR, and EGFR pathways.⁵² These findings suggest that neddylation may facilitate tumor growth and metastasis by modulating the TME through the regulation of numerous signaling pathways in CRC.

To improve the accuracy of risk stratification in CRC patients, a neddylation-related gene signature (NRGS) was developed. This signature demonstrated high predictive accuracy, with patients exhibiting higher risk scores experiencing significantly poorer outcomes across the training, validation, and combined cohorts. The area under the curve (AUC) for 1-year, 2-year, and 3-year overall survival (OS) in the training cohort was 0.979, 0.989, and 0.996, respectively. The NRGS-based risk score was associated with an increased likelihood of recurrence or metastasis, more advanced disease stages, and functioned as an independent prognostic factor for overall survival (OS) in CRC. Genetic mutation analysis indicated that patients with elevated risk scores exhibited a higher mutation frequency compared to those with lower risk scores, suggesting that an increase in genetic mutations may lead to cellular physiological dysfunction and promote tumor metastasis. Notably, the group characterized by elevated NRGS demonstrated a higher prevalence of immunosuppressive cells, including myeloid-derived suppressor cells (MDSCs), mesenchymal stem cells (MSCs), and M2-tumor-associated macrophages (TAMs). Furthermore, this group exhibited an upregulation of chemokines and their corresponding receptors. The intricate interactions among dysregulated cytokines, chemokines, growth factors, and matrix-remodeling enzymes play a critical role in to the pathogenesis of CRC and elicit systemic responses that influence disease outcomes.⁵³ The upregulation of CCR/CCL5 expression has been correlated with the infiltration of immune cells that are associated with poor prognosis, including regulatory T cells (Tregs), M1 and M2 macrophages, myeloid-derived suppressor cells, and cancer-associated fibroblasts.⁵⁴ Additionally, increased CCR/CCL5 expression levels are associated with a wide array of immunosuppressive proteins, such as PD-1, PD-L1, PD-L2, CTLA4, CD80, CD86, TIM3, IDO1, LAG3, and IFN-γ, suggesting potential mechanisms by which CRC circumvents anti-cancer immune responses.⁵⁴ Macrophage-derived CCL5 facilitates the immune evasion of colorectal cancer cells via the p65/STAT3-CSN5-PD-L1 pathway.⁵⁵ Moreover, the neddylation pathway enhances CCL2 transactivation and tumor-associated macrophage (TAM) infiltration,⁵⁶ which are critical for tumor immunosuppression, thereby creating a microenvironment that supports to cancer development and progression.

Significantly, the cohort exhibiting elevated NRGS demonstrated a low IPS score, a low TIDE score, and increased scores for tumor immune exclusion and dysfunction. The IPS has been identified as a superior predictor of response to anti-CTLA4 and anti-PD1 antibodies, with a low IPS score suggesting a diminished response to immunotherapy.⁵⁷ Conversely, a low TIDE score indicates a reduced likelihood of immune escape and an enhanced response to immunotherapy.^58,59 Furthermore, the high NRGS group exhibited resistance to standard colorectal cancer therapies, including oxaliplatin, 5-Fluorouracil, KRAS G12C inhibitors, and sorafenib. Therefore, this signature could help clinicians identify patient subgroups for tailored immunotherapy and chemotherapy, but larger patient groups are needed for validation.

Limitations of the Study

This study represents the inaugural effort to identify prognostically significant neddylation-related genes in colorectal cancer (CRC) through the utilization of scRNA-seq and bulk RNA sequencing data. Nonetheless, several limitations are evident: firstly, the sample size for scRNA-seq was insufficient, and there was notable heterogeneity between scRNA-seq and bulk RNA-seq data. Secondly, although the Neddylation-Related Gene Signature (NRGS) demonstrated efficacy in predicting prognosis in CRC patients, further validation with larger and more diverse clinical samples is necessary to enhance the reliability of the model’s predictive capacity. Additionally, the study’s reliance on public databases without the inclusion of clinical samples, the absence of validation of gene expression in CRC tissues, and the lack of cellular or animal experiments to evaluate the genes’ impact on CRC behavior, constitute significant limitations. Consequently, further research is warranted to elucidate the roles and mechanisms of these genes in CRC progression. Despite these constraints, the identified markers provide valuable guidance for future investigations and hold potential for informing novel therapeutic targets and clinical strategies for CRC.

Conclusion

This study conducted a comprehensive investigation into the roles of neddylation-related genes (NRGs) in the progression of colorectal cancer (CRC), as well as their impact on responses to immunotherapy and chemotherapy. Our results indicated that patients with a high neddylation score exhibited a greater presence of malignant cells and diverse immune cell populations, alongside activation of critical pathways associated with tumor proliferation and immune evasion. We developed a neddylation-related gene signature (NRGS), which effectively characterizes the immunological landscape of CRC patients and provides a reliable prognostic tool. This signature can aid clinicians in identifying specific patient subgroups that may derive benefit from tailored immunotherapy and chemotherapy regimens. Future research should prioritize the validation of these findings in larger patient cohorts and investigate the therapeutic potential of targeting NRGs in the treatment of CRC.

Data Sharing Statement

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Ethics Approval and Consent to Participate

All relevant ethical regulations were followed the original study of the datasets and the authors of the source studies had also obtained informed consent from participants. Ethical approval for this study was exempted by the Fourth Affiliated Hospital of Guangzhou Medical University, as the data were obtained from public sources.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This study were supported by Guangzhou Zengcheng District Science and Technology Innovation fund project (ZCKJ2019-008 and 2021ZCMS18), and Science and Technology Program of Guangzhou, China (2024A03J0948).

Disclosure

The authors declare that they have no competing interests in this work.

References

1. Siegel RL, Wagle NS, Cercek A, Smith RA, Jemal A. Colorectal cancer statistics, 2023. CA: a Cancer Journal for Clinicians. 2023;73(3):233–254. doi:10.3322/caac.21772

2. Ganesh K, Stadler ZK, Cercek A, et al. Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol. 2019;16(6):361–375. doi:10.1038/s41575-019-0126-x

3. Ganesh K. Optimizing immunotherapy for colorectal cancer. Nature Reviews. Gastroenterology & Hepatology. 2022;19(2):93–94. doi:10.1038/s41575-021-00569-4

4. Zhou C, Cheng X, Tu S. Current status and future perspective of immune checkpoint inhibitors in colorectal cancer. CANCER LETT. 2021;521:119–129. doi:10.1016/j.canlet.2021.07.023

5. Zhang S, Yu Q, Li Z, Zhao Y, Sun Y. Protein neddylation and its role in health and diseases. Signal Transduct Target Ther. 2024;9(1):85. doi:10.1038/s41392-024-01800-9

6. Enchev RI, Schulman BA, Peter M. Protein neddylation: beyond cullin-RING ligases. Nat Rev Mol Cell Biol. 2015;16(1):30–44. doi:10.1038/nrm3919

7. Li L, Wang M, Yu G, et al. Overactivated neddylation pathway as a therapeutic target in lung cancer. J Natl Cancer Inst. 2014;106(6):u83. doi:10.1093/jnci/dju083

8. Jones TM, Carew JS, Bauman JE, Nawrocki ST. Targeting NEDDylation as a Novel Approach to Improve the Treatment of Head and Neck Cancer. Cancers (Basel). 2021;13(13):3250. doi:10.3390/cancers13133250

9. Norton JP, Augert A, Eastwood E, Basom R, Rudin CM, MacPherson D. Protein neddylation as a therapeutic target in pulmonary and extrapulmonary small cell carcinomas. GENE DEV. 2021;35(11–12):870–887. doi:10.1101/gad.348316.121

10. Kuo KL, Ho I-L, Shi C-S. MLN4924, a novel protein neddylation inhibitor, suppresses proliferation and migration of human urothelial carcinoma: in vitro and in vivo studies. CANCER LETT. 2015;363(2):127–136. doi:10.1016/j.canlet.2015.01.015

11. Jin Y, Zhang P, Wang Y. Neddylation Blockade Diminishes Hepatic Metastasis by Dampening Cancer Stem-Like Cells and Angiogenesis in Uveal Melanoma. CLIN CANCER RES. 2018;24(15):3741–3754. doi:10.1158/1078-0432.CCR-17-1703

12. Yao WT, Wu J-F, Yu G-Y. Suppression of tumor angiogenesis by targeting the protein neddylation pathway. CELL DEATH DIS. 2014;5(2):e1059. doi:10.1038/cddis.2014.21

13. Luo Z, Yu G, Lee HW. The Nedd8-activating enzyme inhibitor MLN4924 induces autophagy and apoptosis to suppress liver cancer cell growth. CANCER RES. 2012;72(13):3360–3371. doi:10.1158/0008-5472.CAN-12-0388

14. Yang Z, Zhang J, Lin X. Inhibition of neddylation modification by MLN4924 sensitizes hepatocellular carcinoma cells to sorafenib. ONCOL REP. 2019;41(6):3257–3269. doi:10.3892/or.2019.7098

15. Li H, Zhou W, Li L. Inhibition of Neddylation Modification Sensitizes Pancreatic Cancer Cells to Gemcitabine. NEOPLASIA. 2017;19(6):509–518. doi:10.1016/j.neo.2017.04.003

16. Shapiro DD, Zacharias NM, Tripathi DN. Neddylation inhibition sensitises renal medullary carcinoma tumours to platinum chemotherapy. Clin Transl Med. 2023;13(5):e1267. doi:10.1002/ctm2.1267

17. Wei D, Li H, Yu J. Radiosensitization of human pancreatic cancer cells by MLN4924, an investigational NEDD8-activating enzyme inhibitor. CANCER RES. 2012;72(1):282–293. doi:10.1158/0008-5472.CAN-11-2866

18. Vanderdys V, Allak A, Guessous F. The Neddylation Inhibitor Pevonedistat (MLN4924) Suppresses and Radiosensitizes Head and Neck Squamous Carcinoma Cells and Tumors. MOL CANCER THER. 2018;17(2):368–380. doi:10.1158/1535-7163.MCT-17-0083

19. Zhou L, Jiang Y, Luo Q, Li L, Jia L. Neddylation: a novel modulator of the tumor microenvironment. MOL CANCER. 2019;18(1):77. doi:10.1186/s12943-019-0979-1

20. Mao H, Lin X, Sun Y. Neddylation Regulation of Immune Responses. Research (Washington). 2023;6:283.

21. Zhou S, Zhao X, Yang Z. Neddylation inhibition upregulates PD‐L1 expression and enhances the efficacy of immune checkpoint blockade in glioblastoma. INT J CANCER. 2019;145(3):763–774. doi:10.1002/ijc.32379

22. Cui Y, Chen Z, Pan B. Neddylation pattern indicates tumor microenvironment characterization and predicts prognosis in lung adenocarcinoma. Front Cell Dev Biol. 2022;10:979262. doi:10.3389/fcell.2022.979262

23. Li Y, Niu J, Wang Y. Machine learning-based neddylation landscape indicates different prognosis and immune microenvironment in endometrial cancer. Frontiers in Oncology. 2023;13:1084523. doi:10.3389/fonc.2023.1084523

24. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. NAT METHODS. 2008;5(7):621–628. doi:10.1038/nmeth.1226

25. Tang F, Barbacioru C, Wang Y. mRNA-Seq whole-transcriptome analysis of a single cell. NAT METHODS. 2009;6(5):377–382. doi:10.1038/nmeth.1315

26. Tang Q, Tang L, Wang X. Comprehensive Analyses of Single-Cell and Bulk RNA Sequencing Data From M2 Macrophages to Elucidate the Immune Prognostic Signature in Patients with Gastric Cancer Peritoneal Metastasis. Immunotargets Ther. 2025;14:383–402. doi:10.2147/ITT.S506143

27. Tao Y, Li J, Pan J. Integration of scRNA-Seq and Bulk RNA-Seq Identifies Circadian Rhythm Disruption-Related Genes Associated with Prognosis and Drug Resistance in Colorectal Cancer Patients. Immunotargets Ther. 2025;14:475–489. doi:10.2147/ITT.S499806

28. Asada K, Takasawa K, Machino H. Single-Cell Analysis Using Machine Learning Techniques and Its Application to Medical Research. Biomedicines. 2021;9(11):1513. doi:10.3390/biomedicines9111513

29. Sun Y, Wu L, Zhong Y. Single-cell landscape of the ecosystem in early-relapse hepatocellular carcinoma. CELL. 2021;184(2):404–421. doi:10.1016/j.cell.2020.11.041

30. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. BIOINFORMATICS. 2010;26(12):1572–1573. doi:10.1093/bioinformatics/btq170

31. Liu Z, Liu L, Weng S. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. NAT COMMUN. 2022;13(1):816. doi:10.1038/s41467-022-28421-6

32. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC BIOINFORMATICS. 2013;14(7). doi:10.1186/1471-2105-14-7

33. Mayakonda A, Lin D, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. GENOME RES. 2018;28(11):1747–1756. doi:10.1101/gr.239244.118

34. Wu T, Hu E, Xu S. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3):100141. doi:10.1016/j.xinn.2021.100141

35. Zeng D, Ye Z, Shen R. IOBR: multi-Omics Immuno-Oncology Biological Research to Decode Tumor Microenvironment and Signatures. Front Immunol. 2021;12:687975. doi:10.3389/fimmu.2021.687975

36. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform. 2021;22(6). doi:10.1093/bib/bbab260

37. Yang C, Guo Y, Wu Z, Huang J, Xiang B. Comprehensive Analysis of Cuproptosis-Related Genes in Prognosis and Immune Infiltration of Hepatocellular Carcinoma Based on Bulk and Single-Cell RNA Sequencing Data. Cancers (Basel). 2022;14(22):5713. doi:10.3390/cancers14225713

38. Olaizola P, Lee-Law PY, Fernandez-Barrena MG. Targeting NAE1-mediated protein hyper-NEDDylation halts cholangiocarcinogenesis and impacts on tumor-stroma crosstalk in experimental models. J HEPATOL. 2022;77(1):177–190. doi:10.1016/j.jhep.2022.02.007

39. Li K, Shi H, Zhang B. Myeloid-derived suppressor cells as immunosuppressive regulators and therapeutic targets in cancer. Signal Transduction and Targeted Therapy. 2021;6(1):362. doi:10.1038/s41392-021-00670-9

40. Tanaka A, Sakaguchi S. Regulatory T cells in cancer immunotherapy. CELL RES. 2017;27(1):109–118. doi:10.1038/cr.2016.151

41. Seino K, Taniguchi M. Functionally distinct NKT cell subsets and subtypes. The Journal of Experimental Medicine. 2005;202(12):1623–1626. doi:10.1084/jem.20051600

42. Hong SM, Lee A-Y, Kim B-J. NAMPT-Driven M2 Polarization of Tumor-Associated Macrophages Leads to an Immunosuppressive Microenvironment in Colorectal Cancer. Adv Sci (Weinh). 2024;11(14):e2303177. doi:10.1002/advs.202303177

43. Huang C, Ou R, Chen X. Tumor cell-derived SPON2 promotes M2-polarized tumor-associated macrophage infiltration and cancer progression by activating PYK2 in CRC. J Exp Clin Cancer Res. 2021;40(1):304. doi:10.1186/s13046-021-02108-0

44. Han S, Wang W, Wang S. Tumor microenvironment remodeling and tumor therapy based on M2-like tumor associated macrophage-targeting nano-complexes. THERANOSTICS. 2021;11(6):2892–2916. doi:10.7150/thno.50928

45. Liu C, Chikina M, Deshpande R. Treg Cells Promote the SREBP1-Dependent Metabolic Fitness of Tumor-Promoting Macrophages via Repression of CD8(+) T Cell-Derived Interferon-gamma. IMMUNITY. 2019;51(2):381–397. doi:10.1016/j.immuni.2019.06.017

46. Prossomariti A, Piazzi G, Alquati C, Ricciardiello L. Are Wnt/beta-Catenin and PI3K/AKT/mTORC1 Distinct Pathways in Colorectal Cancer? Cell Mol Gastroenterol Hepatol. 2020;10(3):491–506. doi:10.1016/j.jcmgh.2020.04.007

47. Calon A, Espinet E, Palomo-Ponce S. Dependency of colorectal cancer on a TGF-beta-driven program in stromal cells for metastasis initiation. CANCER CELL. 2012;22(5):571–584. doi:10.1016/j.ccr.2012.08.013

48. Jackstadt R, van Hooff SR, Leach JD. Epithelial NOTCH Signaling Rewires the Tumor Microenvironment of Colorectal Cancer to Drive Poor-Prognosis Subtypes and Metastasis. CANCER CELL. 2019;36(3):319–336. doi:10.1016/j.ccell.2019.08.003

49. Pennel K, Hatthakarnkul P, Wood CS. JAK/STAT3 represents a therapeutic target for colorectal cancer patients with stromal-rich tumors. J Exp Clin Cancer Res. 2024;43(1):64. doi:10.1186/s13046-024-02958-4

50. Kotsiliti E. MAPK inhibition in BRAFV600E CRC. Nature Reviews. Gastroenterology & Hepatology. 2023;20(4):201. doi:10.1038/s41575-023-00756-5

51. Zhu G, Pei L, Xia H, Tang Q, Bi F. Role of oncogenic KRAS in the prognosis, diagnosis and treatment of colorectal cancer. MOL CANCER. 2021;20(1):143. doi:10.1186/s12943-021-01441-4

52. Liu D, Che X, Wu G. Deciphering the role of neddylation in tumor microenvironment modulation: common outcome of multiple signaling pathways. Biomark Res. 2024;12(5). doi:10.1186/s40364-023-00545-x

53. Bhat AA, Nisar S, Singh M. Cytokine- and chemokine-induced inflammatory colorectal tumor microenvironment: emerging avenue for targeted therapy. Cancer Commun (Lond). 2022;42(8):689–715. doi:10.1002/cac2.12295

54. Schlechter BL, Stebbing J. CCR5 and CCL5 in metastatic colorectal cancer. J IMMUNOTHER CANCER. 2024;12(5):e008722. doi:10.1136/jitc-2023-008722

55. Liu C, Yao Z, Wang J. Macrophage-derived CCL5 facilitates immune escape of colorectal cancer cells via the p65/STAT3-CSN5-PD-L1 pathway. CELL DEATH DIFFER. 2020;27(6):1765–1781. doi:10.1038/s41418-019-0460-0

56. Zhou L, Jiang Y, Liu X. Promotion of tumor-associated macrophages infiltration by elevated neddylation pathway via NF-kappaB-CCL2 signaling in lung cancer. ONCOGENE. 2019;38(29):5792–5804. doi:10.1038/s41388-019-0840-4

57. Charoentong P, Finotello F, Angelova M. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. CELL REP. 2017;18(1):248–262. doi:10.1016/j.celrep.2016.12.019

58. Fu J, Li K, Zhang W. Large-scale public data reuse to model immunotherapy response and resistance. GENOME MED. 2020;12(1):21. doi:10.1186/s13073-020-0721-z

59. Jiang P, Gu S, Pan D. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. NAT MED. 2018;24(10):1550–1558. doi:10.1038/s41591-018-0136-1