Among the 432,231 admissions recorded in the MIMIC-IV database, 73,181 involved ICU admissions, with 50,920 being first-time ICU admissions. This study focused on patients aged 18 and older diagnosed with SAR based on the MIMIC-IV (v2.2) database. According to the “Rhabdomyolysis: An American Association for the Surgery of Trauma Critical Care Committee Clinical Consensus Document” (23), rhabdomyolysis is defined as creatine kinase (CK) levels exceeding five times the upper limit of normal or over 1,000 U/L. Due to CK typically reaches its peak 24–72 h after muscle injury (24), this study enrolled patients diagnosed with sepsis who developed rhabdomyolysis within 3 days of diagnosis. Sepsis is defined according to the Sepsis-3 criteria or an increased Sequential Organ Failure Assessment (SOFA) score of ≥2. The detailed exclusion criteria are outlined in Figure 1.

Exposure to β-blockers was defined as a prescription containing a β-blocker during the ICU stay. The β-blockers included “Acebutolol,” “Atenolol,” “Esmolol,” “Betaxolol,” “Bisoprolol,” “Labetalol,” “Metoprolol,” “Nadolol,” and “Propranolol.” The routes of administration of β-blocker include intravenous push and oral administration.

Eligible patients were classified into two groups according to β-blocker exposure: those who received β-blocker use after ICU admission (β-blocker group) and those who did not receive β-blockers during their ICU stay (non-β-blocker group).

The covariables were selected based on clinical a priori knowledge as well as their statistical relevance. Covariates with initial regression coefficient changes exceeding 10% were selected by stepwise adjustments in the basic and full models. Collinearity was assessed using the variance inflation factor (VIF), with VIF > 2 indicating collinearity (25, 26). Details on covariates and missing data rates are provided in Supplementary Table 1 and Supplementary Figure 1. Missing data were addressed using multiple imputations via chained equations prior to propensity score matching. All variables presented in Supplementary Table 2 were incorporated into the PSM procedure.

The primary outcome assessed was all-cause in-hospital mortality, while secondary outcomes included 28-day mortality, 90-day mortality, and ICU mortality.

Baseline characteristics are presented as mean ± standard deviation (SD) or median and interquartile range (IQR) for continuous variables, and as number (percentage) for categorical variables. The normality of continuous variables was assessed using Student's t-test or Wilcoxon rank-sum test. Between-group differences in categorical variables were evaluated using Pearson's chi-square test or Fisher's exact test.

To explore the association between β-blocker use and in-hospital mortality, we employed PSM to adjust for potential confounders (Supplementary Figures 2, 3). Variables included in the PSM were selected based on a combination of clinical a priori knowledge and their statistical relevance to both exposure and outcome. To ensure comprehensive adjustment for baseline differences, all baseline covariates presented in Supplementary Table 2 were incorporated into the propensity score estimation. This approach allowed us to account for clinically meaningful confounders while minimizing residual confounding.

For the PSM analysis, propensity scores (PS) were estimated using a logistic regression model that incorporated all baseline variables listed in Supplementary Table 2. These variables were selected based on clinical a priori relevance and their potential statistical association with both treatment assignment and outcomes. Patients in the β-blocker and non-β-blocker groups were matched 1:1 using nearest neighbor matching with a caliper width of 0.1 SD or less. Covariate balance before and after matching was assessed using the absolute standardized mean difference (SMD), with an SMD > 0.1 indicating residual imbalance.

To rigorously estimate the association between β-blocker use and clinical outcomes while minimizing confounding bias, we employed a series of analytical models based on propensity scores (PS). Full Cohort Analyses: Using the entire sample, we performed (a) a PS-adjusted model (including PS as a covariate), (b) Logistic regression in the 1:1 PS-matched cohort, (c) an inverse probability of treatment weighting (IPTW) analysis to balance baseline covariates. The PS was estimated using a multivariable logistic regression model. Matched cohort: Within matched cohort, we performed, (d) a crude logistic regression, and (e) a multivariable-adjusted model for residual confounding. All results are reported as odds ratios (ORs) with 95% confidence intervals (CIs) and are presented together in Table 1.

We conducted subgroup analyses to assess the association between β-blocker use and all-cause mortality within subgroups defined by gender, race, age (<65 years or ≥65 years), Sofa, age and sofa, continuous renal replacement therapy (CRRT), mechanical ventilation support, Diabetic, Renal disease, and Chronic pulmonary disease.

To further validate the association between the use of β-blocker and all-cause mortality, we conducted a sensitivity analysis after excluding patients with myocardial infarction, heart failure, cerebrovascular disease, and peripheral vascular disease. Regarding other potential etiologies of rhabdomyolysis, we additionally performed a sensitivity analysis excluding patients with trauma, intoxication, or seizures, which are common non-infectious causes of rhabdomyolysis. To mitigate potential immortal-time bias, we conducted a series of analyses in the full cohort. Landmark analyses were conducted at 24 and 48 h following ICU admission. Patients were stratified according to whether β-blocker use was initiated before or after each predefined landmark. In addition, a time-dependent stratified Cox proportional hazards model was employed to examine the association between β-blocker use and in-hospital mortality. Patients were categorized into three temporal strata based on the timing of β-blocker initiation: within 24 h, between 24–48 h, and beyond 48 h. β-blocker exposure was modeled as a time-dependent covariate, with separate coefficients estimated for each stratum. Furthermore, sensitivity analyses were performed using multivariable Cox regression among patients who received β-blockers within 24 h and those who initiated use within 48 h.

Data analysis was performed using packages R4.1.2 (The R Foundation1) and Free Statistics software version 1.9. Statistical significance was defined as p < 0.05.

Among the 50,920 patients admitted to the ICU for the first time, a total of 1,194 patients were included, with 563 in the non-β-blocker group and 631 in the β-blocker group before matching. After 1:1 propensity score matching, 292 pairs were successfully matched (Supplementary Table 2). Compared to the non-β-blockers group, patients in the β-blockers group were older (64.63 [15.55] years vs. 54.33 [18.19] years), had a higher mean (SD) Charlson comorbidity index score (5.56 [2.62] vs. 3.98 [3.04]), and had longer mechanical ventilation duration (2.16 [2.99] vs. 1.99 [2.93]). However, they had a lower mean sequential organ failure assessment score (6.81 [4.19] vs. 8.2 [4.72]), lower VIS scores (8.04 [23.94] vs. 12.61 [23.83]), lower CRRT) utilization rate (60 [9.5] vs. 83 [14.7]), and lower mechanical ventilation utilization rate (331 [52.5] vs. 359 [63.8]). After 1:1 propensity score matching, the matched cohort included 584 patients, with comparability between the two groups on all specific variables where SMD was less than 0.1 (e.g., mean [SD] age, 59.68 [16.13] years vs. 59.58 [17.39] years; SMD = 0.006). In the overall matched cohort, MICU patients accounted for the largest proportion of ICU admissions (215 [36.8%]); 352 patients (60.27%) required mechanical ventilation support upon ICU admission. Before PSM, there were notable differences in the distribution of infection sites between the β-blocker group and the non-β-blocker group. The standardized mean differences (SMDs) ranged from 0.009 to 0.1, with bacteremia (2.1% vs. 3.7%; SMD = 0.1), pneumonia (30.1% vs. 27.4%; SMD = 0.061), and urinary tract infection (17.9% vs. 15.5%; SMD = 0.066) showing the largest imbalances. After 1:1 propensity score matching, the distributions of infection sites were well-balanced between the matched cohorts, with all SMDs were below 0.1. Specifically, the SMDs for bacteremia (3.4% vs. 4.1%; SMD = 0.036), pneumonia (32.9% vs. 30.5%; SMD = 0.052), and urinary tract infection (19.5% vs. 18.8%; SMD = 0.017) indicated excellent covariate balance after matching. As shown in Supplementary Table 3, before matching, the β-blocker group showed significantly lower in-hospital mortality compared to the non-β-blocker group (15.4% vs. 30.6%, p < 0.001). After matching, the β-blocker group continued to demonstrate significantly lower in-hospital mortality (16.8% vs. 26.4%, p = 0.009). Similar reductions were observed for 28-day, 90-day mortality and ICU mortality. These findings suggest that β-blocker use is associated with lower mortality in SAR patients even after balancing baseline variables.

In the entire cohort, metoprolol was used significantly more frequently in survivors than in non-survivors (54.1% vs. 33.5%, p < 0.001), while atenolol and betaxolol were used exclusively in survivors. In contrast, esmolol was used more frequently in non-survivors (4.5% vs. 1.7%, p = 0.009). These trends persisted in the matched cohort (metoprolol: 48.7% vs. 33.3%, p = 0.002; esmolol: 6.3% vs. 1.3%, p = 0.004). The time to first β-blocker administration was significantly longer in survivors (median: 4.1 vs. 0 h in the entire cohort; 4.4 vs. 0 h in the matched cohort; both p < 0.001), and the total duration of medication was also significantly longer (median: 48.0 vs. 0 h in the entire cohort; 38.0 vs. 0 h in the matched cohort; both p < 0.001). Regarding the route of administration, oral/nasogastric administration was more common in survivors (28.9% vs. 5.9% in the entire cohort), while intravenous-only administration was more frequent in non-survivors (19.0% vs. 4.5%, p < 0.001). In terms of infection sites, bacteremia was observed only in survivors in the matched cohort (5.0% vs. 0%, p = 0.007), and urinary tract infections were also significantly more common in survivors (21.4% vs. 11.1%, p = 0.009; Supplementary Table 4).

In the propensity score–adjusted model, β-blocker use was associated with a significantly lower odds of mortality (OR = 0.49, 95% CI 0.35–0.69; p < 0.001). Similar findings were observed in the propensity score–matched cohort, in which the odds of death remained significantly lower among β-blocker users (OR = 0.56, 95% CI 0.38–0.84; p = 0.005). The inverse probability of treatment weighting (IPTW) analysis yielded nearly identical results (OR = 0.56, 95% CI 0.43–0.72; p < 0.001). The consistency of effect estimates across the adjusted, matched, and weighted models supports the robustness of the association between β-blocker use and lower mortality.

In the matched cohort, the in-hospital mortality rate was significantly lower in the β-blockers group compared to the non-β-blockers group (49 [16.8%] vs. 77 patients [26.4%]; p = 0.005; Table 1, Figure 2). However, patients in the β-blockers group had a longer median (IQR) length of stay (LOS) in the ICU (4.0 [2.2, 7.8] days vs. 3.8 [2.1, 7.1] days; p = 0.007) and in the hospital (9.1 [5.7, 15.1] days vs. 7.5 [4.8, 12.7] days; p = 0.002; Supplementary Table 1). Univariate logistic regression analysis indicated that the length of ICU stay (LOS ICU) was not significantly associated with mortality (OR, 1.02; 95% CI, 0.99–1.05; p = 0.193). In contrast, the length of hospital stay (LOS hospital) was significantly associated with a lower mortality risk (OR, 0.91; 95% CI, 0.88–0.95; p < 0.001; Supplementary Table 5). Additionally, ICU mortality, 28-day mortality, and 90-day mortality were all significantly lower in the β-blockers group than in the non-β-blockers group (Table 1, Supplementary Figure 4). In the entire cohort and matched cohort, multivariable Cox regression analysis revealed that patients receiving β-blockers had significantly lower rates of in-hospital mortality, ICU mortality, as well as 28-day and 90-day mortality, compared to those not receiving β-blockers (Table 2, Supplementary Figure 5). Given the potential alternative etiologies of rhabdomyolysis, we conducted a Cox multivariable regression analysis after excluding patients with trauma, intoxication, or seizure to demonstrate the robustness of the findings (Supplementary Table 6). Additionally, a sensitivity analysis, which excluded patients with heart failure and myocardial infarction, confirmed that β-blocker use was associated with lower mortality in patients with SAR (Supplementary Table 7).

In patients with SAR, we performed subgroup analyses to assess the interaction between β-blockers and all-cause mortality, as demonstrated in Figure 3 and Supplementary Table 8. In the entire cohort, the adjusted HR for the primary outcome was 0.31 (95% CI: 0.24–0.40). Significant interactions were observed for SOFA score (p < 0.001) and the combined Age–SOFA stratification (p = 0.003), but not for age alone (p = 0.465). Compared with the subgroup with SOFA <6 (HR: 0.78, 95% CI: 0.33–1.80), the use of β-blockers was associated with lower mortality in the subgroup with SOFA ≥ 6 (HR: 0.24, 95% CI: 0.18–0.33). This pattern was consistent in the combined subgroup, with a significant association observed in patients aged ≥65 years with SOFA scores ≥6 (HR: 0.22, 95% CI: 0.14–0.34). In the matched cohort, the overall effect was attenuated but remained significant (HR: 0.51, 95% CI: 0.36–0.75). Interaction effects for SOFA (p = 0.002) and Age–SOFA (p = 0.006) persisted, while age alone remained non-significant (p = 0.479). In the subgroup with a SOFA score > 6, the use of β-blockers was associated with a lower mortality risk (HR: 0.37, 95% CI: 0.24–0.57). A similar pattern was observed in the combined subgroup analyses, with a significant association identified among patients aged ≥65 years with a SOFA score ≥6 (HR: 0.37, 95% CI: 0.20–0.69). In contrast, no such association was observed in the subgroup with a SOFA score <6 or in the corresponding combined subgroup, and the findings did not reach statistical significance. Given the limited sample size of our cohort, the possibility of bias cannot be excluded. Therefore, larger-scale prospective studies are warranted to validate these observations.

In landmark analyses at 24 and 48 h, β-blocker use was consistently associated with lower in-hospital mortality (Supplementary Table 9 and Supplementary Figure 6). Before the 24-h landmark, in-hospital mortality rates were 3.5% without β-blocker vs. 1.5% with β-blocker (OR = 0.407, 95% CI 0.188–0.881). At 48 h, the corresponding rates were 12.2% and 3.2% (OR = 0.252, 95% CI 0.155–0.409). Among patients surviving beyond each landmark, β-blocker use remained associated with reduced risk, with ORs of 0.478 (95% CI 0.370–0.617) at 24 h and 0.605 (95% CI 0.453–0.807) at 48 h. In time-dependent stratified Cox models, β-blocker use was evaluated across three temporal strata with adjustment for covariates. β-blocker use conferred significant time-dependent benefit: hazard was reduced by 59% in the early phase (HR, 0.41; p = 0.026), by 85% in the intermediate phase (HR, 0.15; p < 0.001), and by 61% in the late phase (HR, 0.39; p < 0.001; Supplementary Table 10 and Supplementary Figure 7). In the multivariable Cox regression restricted to patients who initiated β-blocker use within 24 h, the adjusted hazard ratios were 0.52 for in-hospital mortality, 0.50 for ICU mortality, 0.49 for 28-day mortality, and 0.52 for 90-day mortality (all p < 0.001). In the corresponding analysis limited to β-blocker initiation within 48 h, the association with lower mortality remained consistent, with adjusted hazard ratios ranging from 0.45 to 0.49 across all endpoints (all p < 0.001; Supplementary Table 11).