Baseline characteristics

The information about the demographic and clinical characteristics of a total of 37,606 patients with confirmed COVID-19 were extracted from the electronic medical records of ten hospitals in Henan Province (Fig. 1). Based on the rigorous inclusion and exclusion criteria, 6943 patients who received azvudine and 1202 patients who received Paxlovid were included for further analysis in this study. We subsequently applied 2:1 propensity score matching (PSM) to control for confounders, resulting in the inclusion of 2404 patients in the azvudine group and 1202 patients in the Paxlovid group.

Fig. 1

Flowchart of the Henan cohort study design

The baseline clinical characteristics of the included patients are shown in Table 1. Before matching, compared with azvudine recipients, Paxlovid recipients were more likely to be male; presented with more severe illness on admission; had fewer concomitant antibiotics and systemic steroids; had longer time from diagnosis to treatment; had fewer kidney diseases; had more liver diseases, cardio-cerebral diseases, and autoimmune diseases; had higher glomerular filtration rate and alkaline phosphatase levels; and had lower neutrophil, prothrombin time, cholesterol, albumin, and total bilirubin levels. After matching, the characteristics of the patients were balanced between azvudine and Paxlovid groups (all p > 0.05 and all standardized mean differences (SMD) <0.1) (Supplementary Fig. 1).

Table 1 Baseline characteristics of patients with COVID-19 in the Henan cohort before and after propensity score matchingEfficacy of azvudine versus Paxlovid on all-cause death and composite disease progression

In the Henan cohort, there were 469 all-cause deaths, including 288 in the azvudine group and 181 in the Paxlovid group. Compared with Paxlovid treatment, azvudine therapy was related to an obviously lower risk of all-cause death, as determined by the log-rank test (p = 0.038) (Fig. 2a). After multivariable adjustment using Cox regression analysis, the hazard ratio (HR) for all-cause death in the azvudine group was 0.82 (95% confidence intervals (CI): 0.676–0.987, p = 0.036) versus that in the Paxlovid group (Fig. 2c).

Fig. 2

Kaplan–Meier curves and multivariate Cox proportional hazards regression analysis of patients receiving azvudine treatment versus Paxlovid treatment in the Henan cohort. Cumulative hazard of all‐cause death (a) and composite disease progression (b) assessed by Kaplan–Meier curves. c Hazard ratio of all-cause death and composite disease progression after adjusting for all baseline covariates in Table 1. HR hazard ratio, 95% CI 95% confidence interval, PDs Person-days, Incidence events/per 1000 PDs

The secondary outcome was composite disease progression. A total of 681 secondary outcome events occurred, including 446 in the azvudine group and 235 in the Paxlovid group. The cumulative hazard of composite disease progression was not significantly different between the azvudine and Paxlovid treatments (p = 0.95) (Fig. 2b). Cox analysis revealed that after adjusting the confounding factors, the HR was 1.15 (95% CI: 0.975-1.345, p = 0.097) in the azvudine group compared with the Paxlovid group (Fig. 2c).

Sensitivity analyses and cross-region cohort validation

The dependability of clinical research results could be influenced by various factors, including data quality, the analysis population, and statistical methods. After the initial results were obtained, we speculated that different approaches to handling missing data, varying matching models, or different analyzing populations might have affected our findings. Therefore, we conducted the following three sensitivity analyses for the Henan cohort to demonstrate the robustness of our conclusions.

First, when missing values were addressed via mean imputation (Supplementary Table 1), Kaplan–Meier analysis revealed there was no notable difference in the risk of all-cause death between azvudine and Paxlovid groups (p = 0.053) (Supplementary Fig. 2a). However, Cox regression analysis indicated that patients receiving azvudine had a 21% lower risk of all-cause death than did those with Paxlovid treatment (95% CI: 0.658–0.959, p = 0.016) (Supplementary Fig. 2c). The risk of composite disease progression was consistent with results obtained from the original dataset (Kaplan–Meier analysis: p = 0.78; Cox analysis: HR: 1.11, 95% CI: 0.944–1.299, p = 0.209) (Supplementary Fig. 2b, c).

Second, when a Probit model was used to perform a 1:2 greedy match (Supplementary Table 2), Kaplan–Meier analysis (p = 0.0041) (Supplementary Fig. 3a) and Cox regression analysis (HR: 0.73, 95% CI: 0.603–0.884, p = 0.001) (Supplementary Fig. 3c) indicated that azvudine group was related to a lower risk of all-cause death than that in the Paxlovid group. The risk of composite disease progression was not obviously different between the two groups (Kaplan–Meier analysis: p = 0.91; Cox analysis: HR: 1.11, 95% CI: 0.947–1.305, p = 0.194) (Supplementary Fig. 3b, c).

Third, the results remained robust when the analysis was repeated after patients discharged from the hospital on the first day after admission were excluded (Supplementary Table 3). In this analysis, azvudine was related to the decreased risk of all-cause death versus Paxlovid, as demonstrated by Kaplan–Meier analysis (p = 0.031) (Supplementary Fig. 4a) and Cox analysis (HR: 0.80, 95% CI: 0.657–0.964, p = 0.02) (Supplementary Fig. 4c). Furthermore, no obvious difference in cumulative incidence of composite disease progression was observed between two groups (Kaplan–Meier analysis: p = 0.8; Cox analysis: HR: 1.13, 95% CI: 0.960–1.335, p = 0.139) (Supplementary Fig. 4b, c).

In addition to sensitivity analysis, to demonstrate the generalizability of our results, we collected 3270 hospitalized COVID-19 patients from outside the Henan region of Xinjiang Province. After exclusion and PSM, a total of 79 azvudine recipients and 78 Paxlovid recipients were enrolled (Supplementary Fig. 5) (Supplementary Table 4). Kaplan–Meier analysis suggested no notable difference in the risks of all-cause mortality (p = 0.39) or composite disease progression (p = 0.27) between the two groups (Supplementary Fig. 6a, b). Cox regression analysis indicated that the azvudine group had a lower risk of all-cause death than the Paxlovid group did (HR: 0.53, 95% CI: 0.283–0.989, p = 0.046), and with no obvious difference in the risk of composite disease progression (HR: 0.57, 95% CI: 0.322–1.021, p = 0.059) (Supplementary Fig. 6c).

Safety

We collected data on adverse events (AEs) during follow-up for both the azvudine and Paxlovid groups (Supplementary Table 5). Compared with azvudine recipients, patients in the Paxlovid group had a greater risk of Grade 1 AEs, including increased alanine aminotransferase (ALT) (p = 0.013), hypercholesterolemia (p < 0.001), and increased aspartate aminotransferase (AST) (p = 0.047). With respect to Grade 2 AEs, Paxlovid administration was related to greater risks of increased decreased platelet (p = 0.009), increased creatinine (p = 0.018), and ALT (p = 0.036) than azvudine. For Grade 3 and greater SEs, Paxlovid treatment was related to a higher incidence of decreased lymphocyte count (p < 0.001).

Survival prediction model

The LASSO regression analysis included thirty-eight variables measured at hospital admission (Fig. 3a). The optimal λ value, λ.1 se, was chosen on the basis of fivefold cross-validation (Fig. 3b). Eleven features were identified and used to establish a nomogram for predicting 10-, 20-, and 30-day survival probabilities in patients receiving azvudine (Fig. 3c). In the training set, the area under the receiver operating characteristic curve (AUCs) of our model were 80.5%, 76.4%, and 75.6% for 10-, 20-, and 30-day survival, respectively (Fig. 3d); in the test set, they were 81.2%, 73.9%, and 76%, respectively (Fig. 3e); and in the external validation set, they were 78.3%, 70.3%, and 90.4%, respectively (Supplementary Fig. 7a). For patients receiving Paxlovid, an additional nomogram was constructed to predict the 10-, 20-, and 30-day survival probabilities (Supplementary Fig. 8a–c). In the training set, the AUCs of the model were 80.7%%, 83.2%, and 86.1% for 10-, 20-, and 30-day survival, respectively (Supplementary Fig. 8d); in the test set, they were 80.1%%, 81.1%, and 87.8%, respectively (Supplementary Fig. 8e); and in the external validation set, they were 65.7%, 51.9%, and 68.2% (Supplementary Fig. 7b).

Fig. 3

Nomogram for survival prediction of hospitalized COVID-19 patients receiving azvudine in the Henan cohort. a LASSO coefficient profile of 40 features. b Tuning parameter (λ) selection in the LASSO model via fivefold cross-validation. c Nomogram to estimate 10-, 20-, and 30-day survival for COVID-19 patients receiving azvudine. The prediction performance of the LASSO model evaluated by ROC curves in the training set (d) and test set (e)

The concordance index (C-index) of survival prediction model for azvudine recipients was 0.821 in the training set, 0.819 in the test set (Supplementary Fig. 9a, b), and 0.833 in the cross-region validation cohort (Supplementary Fig. 7c). In addition, the C-index for the Paxlovid nomogram was 0.793 in the training set, 0.792 in the test set (Supplementary Fig. 9c, d), and 0.703 in the external validation set (Supplementary Fig. 7d). The calibration curves for both nomograms proved great consistency between the predicted and observed survival probabilities in the training and test sets. The decision curve analysis (DCA) curves indicated that using the nomogram to predict survival in COVID-19 patients receiving azvudine or Paxlovid provided greater net benefit than either predicting all patients as survivors or predicting none as survivors across almost all threshold probabilities in both the training and test sets (Supplementary Fig. 9e–h).

Subgroup analysis

To further investigate the effects of antiviral drugs on clinical outcomes across different populations, we stratified the analysis on the basis of sex, age, severity, vaccination dose, concomitant antibiotic use, systemic steroid use, time from diagnosis to treatment exposure, and comorbidities.

For all-cause death, potentially meaningful interactions suggesting a greater benefit of azvudine over Paxlovid were observed in patients who started treatment >5 days after diagnosis (p for interaction = 0.009, HR = 0.56, 95% CI: 0.39−0.78), in those with primary malignant tumors (p for interaction <0.001, HR: 0.33, 95% CI: 0.20−0.54), and in those without systemic steroid use (p for interaction = 0.004, HR: 0.67, 95% CI: 0.53‒0.84) (Table 2). For composite disease progression, potentially meaningful interactions suggesting a greater benefit of azvudine over Paxlovid were observed in patients with moderate COVID-19 (p for interaction = 0.036, HR: 0.67, 95% CI: 0.45−1.00) and with primary malignant tumors (p for interaction = 0.012, HR: 0.54, 95% CI: 0.33−0.88) (Table 2).

Table 2 Subgroup analyses in the Henan cohort for the effects of different drug treatments on all-cause death and composite outcomes according to baseline characteristicsAntitumor effects of azvudine in vivo and in vitro

Subgroup analysis indicated that azvudine had a more pronounced benefit in reducing clinical outcomes in patients with primary malignant tumors. To further explore the antitumor role of azvudine, we investigated its role in HCC cell lines (H22, Hep3b, and Huh7) and lung cancer (LC) cell lines (NCI-H1975 and NCI-H82). Interestingly, CTG assays revealed that azvudine obviously inhibited the proliferation of all the tested cell lines (Fig. 4a and Supplementary Fig. 10). In contrast, no notable difference was observed in the Paxlovid versus the control. Similarly, colony formation assays demonstrated that azvudine markedly decreased the number of colonies (Fig. 4b). Additionally, invasion assays revealed that, compared with Paxlovid or DMSO treatment, azvudine treatment significantly reduced the number of invading HCC cells (Fig. 4c). These results suggest that azvudine can effectively inhibit cancer cell proliferation and invasion.

Fig. 4

Azvudine inhibits tumor growth in vitro and in vivo. a The proliferation of Huh7, Hep3b, and H22 cells treated with azvudine, Paxlovid, and cisplatin was assessed via the CTG assay. b The proliferation of Hep3b and Huh7 cells treated with azvudine, Paxlovid, or DMSO was assessed by a colony formation assay (n = 3 per group). Scar bar, 10 mm. c The invasion of Hep3b and Huh7 cells treated with azvudine, Paxlovid, or DMSO was evaluated by Transwell assays (n = 3 per group). *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. Scar bar, 200 μm. d Tumor growth and weight were monitored in BALB/c mice subcutaneously injected with H22 cells (n = 6 per group). e H&E and IHC staining of Ki-67 and PCNA expression in tumors (n = 6 per group). The error bars represent the means ± SEs. Scar bar, 200 μm

To further evaluate the antitumor effects of azvudine in vivo, an allograft mouse model was generated by subcutaneously injecting H22 cells into mice. The mice were then randomly assigned into three groups: the azvudine group (1 mg/kg, once daily), the Paxlovid group (61.65 mg/kg nirmatrelvir and 20.55 mg/kg ritonavir, twice daily), and the control group (solvent). In line with the in vitro results, tumor growth rate and average tumor weight in the azvudine group were obviously lower than those in the control, and Paxlovid groups (Fig. 4d). Immunohistochemistry (IHC) analysis revealed azvudine group had a weaker Ki-67 and PCNA staining than in the Paxlovid and control groups (Fig. 4e). With respect to the biosafety of azvudine and Paxlovid, blood biochemical indicators of liver and kidney function in mice were within the normal range across all three groups (Supplementary Fig. 11a). Additionally, H&E staining of the heart, lungs, spleen, liver, and kidneys further confirmed the safety of both agents (Supplementary Fig. 11b). These results validated the antitumor effects of azvudine in HCC.

Single-cell RNA sequencing (scRNA-seq) analysis of the tumor immune microenvironment after azvduine treatment

We subsequently performed scRNA-seq on mouse HCC tumor tissue to explore the potential mechanism of the antitumor effect of azvudine. Nine distinct immune cell types were recognized by uniform manifold approximation and projection (UMAP) analysis (Fig. 5a). The distribution of immune cell types in the control group differed from that in the azvudine group (Fig. 5b). Macrophages, monocytes, and CD4+ T cells composed the dominant immune cell types in the two groups (Fig. 5c). The proportions of CD4+ T cells, CD8+ T cells, B cells, and NK cells were significantly greater in the azvudine group than in the control group (Fig. 5d). Conversely, the proportion of macrophages was decreased in the azvudine group. We analyzed the expression levels of typical molecular markers in both groups to investigate the effect of azvudine on functional alterations in immune cells (Fig. 5e). Following treatment with azvudine, the costimulatory molecules ICOS and TNFRSF9 were upregulated in CD4+ T cells, whereas Tnfrsf11 was downregulated in CD4+ T cells; the cytokine receptor IL17a was upregulated in CD4+ T cells; IL1R1 was upregulated in mast cells; and the inhibitory receptor ENTPD1 was downregulated in granulocytes and monocytes.

Fig. 5

scRNA-seq revealed that the global tumor immune profiles respond to azvudine treatment. a UMAP plot of the cluster analysis of all samples, with each color denoting a distinct cell type. b Distribution of nine immune cell types in the control and azvudine groups (n = 6 per group). c Proportion of immune cell types in the control and azvudine groups (n = 6 per group). d Frequencies of different immune cell types in the control and azvudine groups. *p < 0.05, **p < 0.01, ***p < 0.001. e The expression levels of characteristic genes across nine immune cell types in the control and azvudine groups. f UMAP plot of T lymphocyte subclusters in all samples. g Expression profiles of specific marker genes for T cell subclusters in UMAP. h Proportion of T cell subclusters in the control and azvudine groups. C cluster

To identify the phenotypic heterogeneity of T cells, UMAP analysis were performed on the scRNA-seq data from T cells (Fig. 5f–h). A total of 14T lymphocyte clusters were recognized. Compared with those in the control group, the proportions of subclusters 0 (MT2+CD4+ T cells) and 2 (IKZF2+CD4+ T cells) were notably greater in azvudine group, while those of subclusters 5 (IL18RAP+CD4+ T cells), 6 (LY6C2+CD8+ T cells), and 7 (CXCR6+CD4+ T cells) were significantly lower (Supplementary Fig. 12a). Following azvudine treatment, the effector molecule KLRG1 was downregulated in subcluster 2, GZMC was upregulated in subcluster 13 (Supplementary Fig. 12b), the inhibitory receptor TIGIT was upregulated in subcluster 12, and the costimulatory molecule CCL1 was upregulated in subclusters 4 and 7.

The ligand‒receptor interactions across diverse immune cell types were analyzed via CellChat (Fig. 6a, b). The total interaction strength and intercellular interaction number were lower in the azvudine group than in the controls (Fig. 6c, d and Supplementary Fig. 12c). Ligand‒receptor analysis revealed that azvudine predominantly regulated signal transduction from tumor CD4+ T cells to other immune cells through SPP1-(ITGAV + ITGB5) and SPP1-(ITGAV + ITGB1), as well as COL1A1-(ITGA9 + ITGB1), COL1A1-CD44, and COL1A1-SDC4 (Fig. 6e). Azvudine mainly mediated intercellular interactions from CD8+ T cells to other immune cells through influencing SPP1-(ITGAV + ITGB1), SPP1-(ITGA9 + ITGB1), and SPP1-(ITGA5 + ITGB1). Additionally, azvudine regulated the interactions from other tumor immune cells to CD4+ T and CD8+ T cells through ligand‒receptor pairs such as COL1A2‒(ITGA1 + ITGB1), CCL4‒CCR5, CCL6‒CCR2, and CXCL16‒CXCR6 (Supplementary Fig. 13).

Fig. 6

Azvudine reduces cell–cell interactions and communication among immune cells in HCC. The number of intercellular interactions in the control (a) and azvudine groups (b). Heatmap displaying the numbers of cell–cell interactions in the control (c) and azvudine (d) groups. e Bubble map presenting the comparison of significant ligand–receptor pairs in the control and azvudine groups when CD4+ T and CD8+ T cells were used as ligand cells

The pseudotime trajectories of the CD4+ and CD8+ T cell subclusters between two groups were obtained via the Monocle algorithm (Supplementary Fig. 14a, d). The results revealed that azvudine may influence the development and differentiation of CD4+ and CD8+ T cells (Supplementary Fig. 14b, e). We observed that subclusters 4, 5, 6, and 7 were comprised primarily of cells in the initial stage of development that may differentiate into subclusters 0, 9, 11, 12, and 13 (Supplementary Fig. 14c, f). KEGG enrichment analysis suggested that the differentially expressed genes (DEGs) in CD4+ T cells were involved mainly in cytokine‒cytokine receptor interactions and prion disease (Supplementary Fig. 15a, b), whereas DEGs in CD8+ T cells were involved primarily in protein processing in endoplasmic reticulum and cytokine‒cytokine receptor interactions (Supplementary Fig. 15c, d).