Ongoing immune response and neutrophil changes after COVID-19 with possible treatment ideas

The IAV and SARS-CoV-2 negative P. roborovskii hamsters (6–8 weeks old) were purchased from Dooyeol Biotech. Animal studies were conducted following experimental procedures approved by the Institutional Animal Care and Use Committee of Chungbuk National University (approval number CBNUA-194R-22-01) and the Institute for Basic Science (approval number IBS-2024-001). All experiments were conducted in biosafety level 3 laboratories at Chungbuk National University (KCDC-14-3-07) and the Institute for Basic Science (KCDC-23-3-06).

A mouse-adapted pandemic H1N1 IAV (maCA04/09)¹¹ and a SARS-CoV-2 delta variant strain (GISAID accession number EPI_ISL_10189547) isolated from a patient with confirmed COVID-19 in South Korea were cultured in MDCK (Madin-Darby Canine Kidney cells) and Vero E6 (African green monkey kidney cells) cells, respectively^12,13. For virus infections, groups of hamsters were anesthetized using isoflurane, and administered 30 µl of SARS-CoV-2 at a titer of 10^5.5 tissue culture infective dose 50% (TCID₅₀)/ml or IAV at 10⁷ TCID₅₀/ml through intranasal routes. Body weight and survival were monitored every 2 dpi until the end of experiments. Animals that experienced more than a 25% reduction in body weight were humanely euthanized.

WRW4 (Selleckchem) was used as an FPR2 antagonist, Paquinimod (TargetMol) as an S100A8/A9 inhibitor and Sivelestat (MedChemExpress) as a neutrophil elastase inhibitor. All compounds were dissolved in 2% dimethyl sulfoxide (DMSO) in saline. Drug administration was performed using two schedules: early treatment from 5 to 11 dpi (n = 50) and delayed treatment from 15 to 21 dpi (n = 20) following SARS-CoV-2 infection. With the exception of the negative control and nontreated groups, the following treatments were applied: 2% DMSO in saline (50 μl, intraperitoneally), Paquinimod (10 mg/kg/day, 50 μl, orally) and Sivelestat (10 mg/kg/day, 50 μl, intraperitoneally), administered daily for 7 days. WRW4 (8 mg/kg/day, 50 μl, intraperitoneally) was administered by 2 days for a total of 4 days.

For titration of TCID₅₀, tissue homogenates were prepared by homogenizing entire tissue samples in 0.7 ml of virus culture media containing 2% penicillin–streptomycin (P/S, Gibco) using a TissuLyser II (Qiagen). The homogenates were clarified through three rounds of centrifugation at 12,000 rpm for 10 min. Influenza virus titers were determined by infecting MDCK cells with 10-fold serial dilutions of tissue homogenates or BALF in Eagle's minimum essential medium (Gibco) containing 1% P/S, followed by incubation with Eagle's minimum essential medium supplemented with 1 μg/ml TPCK-trypsin (Worthington Biochemical) at 37 °C with 5% CO₂ for 72 h. Similarly, SARS-CoV-2 titers were assessed by infecting Vero E6 cells with 10-fold serial dilutions of tissue homogenates or BALF in Dulbecco's modified Eagle medium (Gibco) with 1% P/S, followed by incubation with Dulbecco's modified Eagle medium containing 2% fetal bovine serum and 1% P/S at 37°C with 5% CO₂ for 72 h (ref. ¹²). Cytopathic effects were monitored daily, and viral titers were calculated as log₁₀ TCID₅₀/ml using the Reed–Muench method¹⁴.

Harvested hamster tissues were placed in TRIzol (Invitrogen) and homogenized using an Omni Tissue Homogenizer (OMNI International). RNA was extracted from the homogenized tissues following the TRIzol RNA extraction protocol. The extracted RNA was then used to synthesize complementary DNA (cDNA) with SuperScript III Reverse Transcriptase (Invitrogen). For the synthesis, oligo d(T) and random primers (Promega) were used. qRT–PCR was conducted using the synthesized cDNA, IQ SYBR Green SuperMix (Bio-Rad) and primers (Supplementary Table) and was detected by CFX Opus 96 (Bio-Rad). 1

The GM-CSF and G-CSF ELISAs were performed using the Hamster GM-CSF Antibody ELISA Kit (MyBioSource) and G-CSF ELISA Kit (MyBioSource), following the protocols provided by the manufacturer.

Whole-lung samples from mock and virus-infected mice were collected at 5, 15 and 30 dpi. The samples were fixed in 10% neutral-buffered formalin, embedded in paraffin and stained with hematoxylin and eosin (H&E) or Masson's trichrome (MT) staining solution. Immunohistochemistry (IHC) was performed using the Leica Bond Rx Research Stainer (Leica Microsystems) with bond polymer refine detection. For the detection of SARS-CoV-2 antigens, tissue sections were stained with anti-SARS-CoV-2 spike protein S1 antibody (HL134, Abcam) and anti-SARS-CoV-2 nucleocapsid rabbit monoclonal antibody (EPR24334-118, Abcam). Detection was performed using peroxidase AffiniPure goat anti-rabbit IgG (H + L) secondary antibody (Jackson ImmunoResearch). For detection of IAV antigens, sections were stained with an anti-nucleoprotein antibody (H16-L10-4R5, GeneTex), followed by Peroxidase AffiniPure Goat Anti-Mouse IgG (H + L) secondary antibody (Jackson ImmunoResearch).

The RNAscope in situ hybridization assay used V-Influenza-H1N1-H5N1-M, V-nCoV2019-S probes (ACD Bio) and the RNAscope 2.5 HD Detection Kit – RED (ACD Bio), following the manufacturer's protocol. All the stained sections were visualized using an Axio scan Z1 (Carl Zeiss), and fibrosis areas were analyzed with pixel classification in QuPath software (v0.5.0).

The samples processed using the Leica Bond Rx Research stainer (Leica Microsystems) with Opal 3-plex detection kit (PerkinElmer). Anti-myeloperoxidase (MPO) antibody (Abcam) and anti-CD68 antibody (Abcam) were used for detection, followed by peroxidase-conjugated anti-rabbit antibody (Jackson Immuno Research Laboratories). Fluorophore labeling was performed using Opal 570 and 690 fluorophores, and the samples were mounted with Prolong Gold Antifade reagent with DAPI (Invitrogen). Quantitative analysis was conducted using QuPath software (v0.5.0), performing object classification to quantify results.

The single-cell RNA libraries were prepared utilizing the 10x Genomics Chromium Single Cell 3′ Reagent Kit (v3.1, 10x Genomics) based on the manufacturer's instructions. In brief, from each tissue sample, approximately 8,000 cells were loaded in one channel. In each droplet, cDNA was generated through reverse transcription. The quality of the resulting cDNA libraries was evaluated using the 4150 TapeStation 4150 (Agilent) and quantified with the High Sensitivity D1000 Screen Tape (Agilent) based on the manufacturer's library quantification protocol. Following cluster amplification of denatured templates, sequencing was performed in a paired-end fashion (2 × 100 bp) using the Illumina Novaseq 6000 platform (Illumina).

Genomes of SARS-CoV-2 delta variant, H1N1 IAV and P. roborovskii hamsters (NCBI accession number GCF_943737965.1) were retrieved from the NCBI database (as of October 2023). To improve the gene annotation, we used the Eggnog-mapper software (v2.1.11)¹⁵. The scRNA-seq data were processed using Cell Ranger from 10x Genomics (v7.2.0; www.10xgenomics.com↗) and analyzed with the Seurat package (v4.1.1)¹⁶. In brief, following the alignment of sequencing reads¹⁷, low-quality cells were removed on the basis of the unique molecular identifier counts (<400) and gene counts (<200 or >8,000), and the proportion of unique molecular identifier counts mapped to mitochondrial gene/cell (≥20%). In addition, cell-free RNA and doublets were removed using SoupX (v1.6.2) and DoubletFinder (v2.0), respectively^18,19. Qualified scRNA-seq data were then normalized using SCTransform²⁰, reduced in dimensionality by principal component analysis and integrated using reciprocal principal component analysis. The batch-corrected matrices were constructed using the k-nearest neighbor graph with the FindNeighbors function (k = 50). Finally, single-cell transcriptomic networks for three tissue types were visualized using the Uniform Manifold Approximation and Projection (UMAP) method. Clustering was performed using the Louvain algorithm, and cell types were annotated on the basis of canonical differentially expressed genes (DEGs, avg_log₂fold change (FC) >0.25; P < 0.05, Benjamini–Hochberg corrected) and verified with CellKb (v2.6; https://www.cellkb.com/↗). For trajectory analyses of myeloid cell populations, we used the Monocle 3 package (v1.3.5)²¹.

From the scRNA-seq data, the cell type percentages were calculated by dividing the number of cells in each type by the total cell count in the control and diseased groups. To minimize bias from small sample sizes, population changes were quantified by dividing diseased cell counts by controls for all hamster pairs. Following logarithmic formation (log₂), we conducted statistical analyses on the relative variation in cell count between the control and diseased groups using a Wilcoxon rank-sum test.

For the functions of interest in selected cell types, we assessed the combined expressions of a set (module) of relevant genes per cell using the AddModuleScore in the Seurat (v4.1.1)¹⁶. To this end, the Gene Ontology Biological Process (GOBP) database, along with the 50 hallmark gene sets and WikiPathway gene sets from the Molecular Signatures Database (MSigDB, v7.5.1), were used²². Module scores for neutrophil granule production and aging were evaluated on the basis of gene sets informed by prior studies on neutrophil function and aging processes^23,24. We used the Wilcoxon rank-sum test to statistically compare variations in all these module scores between diseased versus control groups, respectively.

Cell-to-cell communication networks were performed using CellChat R package (v2.1.2)²⁵. The lung scRNA-seq dataset was input into the createCellChat function, followed by preprocessing with the identifyOverExpressedGenes and identifyOverExpressedInteractions functions. Communication probabilities were computed using the computeCommunProb function with type = truncatedMean and trim = 0.2 and inferred at signaling pathway levels by the computeCommunProbPathway function. The TGF-β signaling pathway was visualized by the netVisual_aggregate function.

To identify sets of genes that are specifically upregulated or downregulated in the diseased groups, we calculated the FC in the expression of genes and statistical significance (Wilcoxon rank-sum test) in the respective diseased groups compared with the control counterparts. From these comparisons, the group-specific upregulated genes were considered as those that showed the log₂FC >0.25 and P value <0.05 in each diseased group but not in the other two. Conversely, group-specific downregulated genes were defined as those showing the log₂FC <0.25 and P value <0.05 in the given group but not in the others. Functional analyses using these group-specific upregulated or downregulated genes were performed by EnrichR (v3.2) software with the databases of GOBP, MSigDB Hallmark, Kyoto Encyclopedia of Genes and Genomes (KEGG) and/or Reactome²⁶. Gene pathway analysis was conducted with Ingenuity pathway analysis (IPA; Qiagen).

To assess the composition of immune cell populations, lung and spleen tissues were collected from euthanized animals at 5, 15 and 30 dpi. Tissues were dissociated into single-cell suspensions using the gentleMACS Octo Dissociator (Miltenyi Biotec) and filtered through a 70-μm cell strainer. Approximately 1 × 10⁶ cells per sample were incubated with anti-mouse CD16/CD32 antibody (Mouse Fc Block, BD Biosciences) at 4 °C for 30 min to block Fc receptors. After centrifugation and removal of the supernatant, cells were stained at 4 °C for 30 min with the following surface antibodies: BV421-conjugated anti-human CD11b antibody (clone ICRF44, BD Biosciences), PerCP-Cy5.5-conjugated anti-human CD14 antibody (clone M5E2, BD Biosciences) and FITC-conjugated anti-mouse Ly-6G antibody (clone 1A8, BioLegend). Cells were then washed with Flow Cytometry Staining Buffer (Invitrogen) and fixed in 4% paraformaldehyde.

For intracellular staining, fixed cells were permeabilized using the Cytofix/Cytoperm Plus Fixation/Permeabilization Solution Kit (BD Biosciences) at 4 °C for 30 min. After washing twice with Perm/Wash buffer, cells were incubated with PE-conjugated anti-human, mouse, ferret CD79a antibody (clone HM47, Invitrogen) and Alexa Fluor 647-conjugated anti-human CD3 antibody (clone CD3-12, Bio-Rad) at 4 °C for 30 min. Following two final washes, cells were fixed again in 4% paraformaldehyde and analyzed on a FACSymphony A3 Cell Analyzer (BD Biosciences). Flow cytometric data were processed using FlowJo software (v10.9.0).

Data are shown as mean ± standard deviations (s.d.). The statistical analyses of in vitro and in vivo studies were performed using one-way analysis of variance (ANOVA) (Figs. 1, 3 and 5–7 and Supplementary Figs. 1, 6 and 7). For other analyses, the Wilcoxon rank-sum (Mann–Whitney U) test was used to assess statistical significance (Figs. 2–5 and 8 and Supplementary Figs. 2–5). Statistical analyses were performed using GraphPad Prism 10 (v10.3.1). Detailed statistical methods and results are provided in each figure legend.

To investigate whether the observed differences in recovery dynamics were related to variations in viral titers and replication kinetics, respiratory tissues were collected at 5, 15 and 30 dpi. In the SARS-CoV-2 group, infectious virus was recovered at 5 dpi, with lung titers approximately 10 times higher than those in the nasal turbinates, despite similar viral RNA copy numbers in both tissues (Fig. 1h, i). Similarly, in the IAV group, both infectious virus and viral RNA were detected in nasal turbinates and lung tissues at 5 dpi (Fig. 1j, k). However, by 15 and 30 dpi, no infectious virus or viral RNA was detectable in either group (Supplementary Fig. 1a–h). Viral RNA levels in nasal turbinates and lungs were analyzed at two-day intervals from 3 to 15 dpi, revealing higher viral titers in the lungs, peaking at 3 dpi and declining until undetectable by 13 dpi (Fig. 1l, m and Supplementary Fig. 1a–h). These results indicate that the persistent weight loss in the SARS-CoV-2 group is unrelated to ongoing viral replication.

Histopathological examination revealed that virus-induced lung lesions, including fibrosis, persisted until 30 dpi only in the SARS-CoV-2_non-recovery group (Fig. 1n, o). By contrast, the IAV_recovery and SARS-CoV-2_recovery groups showed significant improvement in lung pathology by 30 dpi, as confirmed by H&E and MT staining. The IAV_recovery group showed fibrosis at 5 dpi, which resolved over time, while fibrosis persisted in the SARS-CoV-2_non-recovery group (Fig. 1n, o). RNA in situ hybridization (RNAscope) confirmed viral RNA presence only until 5 dpi in both groups, consistent with TCID₅₀ and qRT–PCR results (Fig. 1h–k, n and Supplementary Fig. 1a–h). These findings suggest that persistent weight loss and severe lung pathology in the SARS-CoV-2_non-recovery group are driven by sustained lung damage rather than prolonged viral replication, highlighting a distinct pathophysiological feature of long COVID in this animal model²⁷.

Next, we conducted scRNA-seq analysis of 158,116 individual cell transcriptomes from 36 tissue samples, including lung, BALF and spleen (Supplementary Table 2). Unsupervised clustering and differential gene expression analysis identified 15 major cell types, including immune, progenitor, epithelial and endothelial cells (Figs. 2d–h and Supplementary Table 3). Cell-type quantification revealed that neutrophil and macrophage (MQ) compositions in both recovery groups and CTRL were comparable across BALF and lung tissues (neutrophils 3.04–23.08%, MQ 6.12–80.05%) (Fig. 2i, j, Supplementary Fig. 2b, c and Supplementary Table 4). In stark contrast, the SARS2_PASC group exhibited a marked increase in neutrophils (61.88–78.45%) and a decrease in MQ (4.16–16.02%) in both tissues. In addition, myeloid progenitor cells notably accumulated in the spleen of the SARS2_PASC group, whereas they were scarce in the recovery groups (Fig. 2k, Supplementary Fig. 2d and Supplementary Table 4).

Beyond these PASC-specific alterations, recovery groups exhibited distinct immune cell dynamics. For example, the SARS2_rec group displayed increased B cell proportions in BALF and spleen compared with CTRL and other disease groups (Fig. 2i, k). Conversely, the IAV_rec group showed a twofold increase in T cell composition in the lung relative to CTRL, a change absent in the SARS2_PASC and SARS2_rec groups (Fig. 2j). These findings highlight unique transcriptional profiles and cellular changes linked to SARS-CoV-2 and IAV infections, emphasizing the distinct immune dysregulation associated with SARS2_PASC.

To investigate the mechanisms underlying the increased neutrophil levels and decreased MQ levels in the SARS2_PASC group, we conducted a trajectory analysis of myeloid cell populations in the UMAP plots to investigate their compositional differences (Fig. 3e). To this end, myeloid progenitors were classified into granulocyte–monocyte progenitor (GMP), neutrophil progenitor (Neu_prog) and monocyte progenitor (Mono_prog) subtypes using eight marker genes (Supplementary Fig. 3a,b). Our trajectory analysis identified two distinct differentiation pathways: GMP to neutrophils via Neu_prog (path1) and GMP to monocytes/MQs via Mono_prog (path2) (Fig. 3e and Supplementary Fig. 3c). Accordingly, neutrophil-specific transcription factors (TFs) were highly expressed in Neu_prog populations, while monocyte/MQ-specific TFs were more prominent in Mono_prog populations (Supplementary Fig. 3d). Notably, in the SARS2_PASC group, the GMP, Neu_prog and neutrophil populations exhibited increased expression of neutrophil-specific TFs, particularly CEBPE, which is critical for neutrophil maturation, along with other key TFs such as KLF5 and GFI1^32–35 (Fig. 3f and Supplementary Fig. 3e, f). Conversely, TFs associated with monocyte/MQ differentiation, such as IRF8³⁶, were substantially downregulated in the SARS2_PASC group compared with recovery groups, which had higher monocyte and MQ ratios (Fig. 3g). These findings suggest that the skewed neutrophil-to-monocyte/MQ ratios observed in the SARS2_PASC group are driven by disrupted myeloid cell differentiation pathways.

Although both SARS-CoV-2 and IAV infection groups effectively suppressed viral replication by 13 dpi, the persistent neutrophil accumulation and reduced MQ population observed exclusively in the SARS2_PASC group remain unexplained. We hypothesized that, despite the absence of detectable viral RNA, other factors might drive sustained tissue inflammation in this group. Recent studies have suggested that the SARS-CoV-2 spike protein, particularly the S1 subunit, can persist in tissues and contribute to prolonged sequelae, even in the absence of viral RNA^37,38. To test this hypothesis, we investigated the presence of viral antigens, specifically the S1 subunit, which is known to modulate key TFs involved in immune response regulation^39,40. IHC revealed that S1 antigens remained detectable in the SARS2_PASC group up to 30 dpi but were absent in the SARS2_rec group (Fig. 3h, i and Supplementary Fig. 3g). By contrast, SARS-CoV-2 nucleocapsid (N) and IAV nucleoprotein (NP) were detectable only up to 5 dpi and completely cleared by 15 dpi (Fig. 3h, j, k and Supplementary Fig. 3g). In the SARS2_PASC group, the S1 antigens were distributed within interstitial pneumonia lesions, accompanied by inflammatory cell infiltration, indicative of ongoing inflammation (Fig. 3h, l and Supplementary Fig. 3g). Notably, monocytes/MQs were observed to internalize S1 antigens, with neutrophils infiltrating the surrounding areas (Fig. 3l, arrowheads and asterisks). These findings highlight a distinct pattern of chronic inflammation in SARS-CoV-2-infected tissues, driven by the prolonged presence of S1 antigens, a phenomenon not observed in IAV-infected tissues, where NP was cleared much earlier. This finding underscores a critical pathophysiological difference between SARS-CoV-2 and IAV infections. While both viruses achieve similar replication suppression timelines, the residual S1 antigens in SARS-CoV-2-infected tissues appear to sustain inflammatory responses, consistent with previous findings^38,41. Moreover, the prolonged presence of S1 antigens may influence TF expression, skewing immune responses and contributing to the non-recovery phenotype observed in the SARS2_PASC group. These results provide key insights into the mechanisms underlying PASC and suggest that persistent viral antigens, rather than active replication, drive chronic inflammation in this condition.

Given the important roles of myeloid populations in lung fibrosis, to further investigate the characteristics of MQ differentiation from myeloid cells during SARS-CoV-2 infection, we classified MQs into three distinct subpopulations: M1 MQs (M1_MQ), M2 MQs (M2_MQ) and alveolar MQs (Alveolar_MQ) (Fig. 4a and Supplementary Fig. 4e). In respiratory tissues such as BALF and lung, alveolar MQs dominated the MQ population, comprising 75.30–91.30% of total MQs. Conversely, monocytes, M1_MQ and M2_MQ were more prevalent in the spleen, highlighting tissue-specific distribution patterns of MQ subsets (Fig. 4a). Comparative functional analyses revealed substantially heightened inflammatory responses in monocytes and MQs of the SARS2_PASC group compared with recovery (SARS2_rec, IAV_rec) and CTRL groups (Fig. 4b, c and Supplementary Table 5). Notably, monocytes, M2_MQ and alveolar MQs in the lungs of the SARS2_PASC group exhibited upregulation of genes associated with TGF-β, VEGF and inflammation pathways, indicating their active involvement in localized inflammation and fibrosis. By contrast, M1_MQ populations in the spleen displayed more pronounced engagement in these pathways, reflecting tissue-specific functional divergence (Fig. 4c).

Importantly, monocytes, M2_MQ and alveolar MQs in the lungs of the SARS2_PASC group persistently overexpressed key extracellular matrix (ECM) genes, including FN1, VCAN, TIMP1 and TIMP2, as well as TGFB1, a critical profibrotic cytokine that stimulates ECM production^42,46 (Fig. 4d and Supplementary Fig. 4f). These genes are closely linked to fibrosis and tissue remodeling, suggesting that these MQ subpopulations are central to sustained ECM deposition and fibrotic progression. Furthermore, fibrosis-associated pathways were highly upregulated in monocytes and M2_MQ populations of the SARS2_PASC group compared with the CTRL group, emphasizing their role in perpetuating lung fibrosis during PASC. In addition, monocytes and M2_MQs in the lungs and spleen of the SARS2_PASC group showed elevated expression of inflammation-related genes, such S100A8, S100A9 and FPR1 (Fig. 4d, e and Supplementary Fig. 4f, g). This suggests that these subpopulations actively contribute to the inflammatory milieu observed in PASC, further differentiating the SARS2_PASC group from the recovery and control groups.

To explore functional differences among neutrophils, we further categorized them into five subclusters (Neu1–Neu5) based on marker genes associated with neutrophil maturation and function (Fig. 5c, d). Cell distribution analysis revealed that neutrophil numbers decreased in both the IAV_rec and SARS2_rec groups compared with the CTRL group (Fig. 5e). Conversely, all neutrophil subclusters increased markedly in the SARS2_PASC group. Unlike their counterparts in the recovery groups, neutrophils in the SARS2_PASC group (Neu1–Neu5) displayed high expression of genes associated with 'leukocyte migration', 'granule production' and 'inflammatory response', including S100A8, S100A9, FPR1, FPR2, MMP8 and MMP9 (Fig. 5f–h and Supplementary Fig. 5a–c). ELISA confirmed elevated levels of GM-CSF and G-CSF, critical for neutrophil development, survival and function, in the SARS2_PASC group^47,48 (Fig. 5i). These increases corresponded with the upregulation of their receptors, CSF2RA and CSF3R, in all Neu1–Neu5 subsets (Fig. 5j). Collectively, these findings indicate that the coordinated upregulation of receptor-cytokine genes supports sustained neutrophil activation and inflammation in the SARS2_PASC group.

Neutrophil trajectory analysis revealed two differentiation pathways: Neu4 via Neu1 and Neu2 (path1) and Neu5 via Neu1 and Neu3 (path2) (Fig. 5k, l and Supplementary Fig. 5d). Neu4, expanding across all tissues, exhibited mature neutrophil characteristics⁴⁹ (Fig. 5d and Supplementary Fig. 5e). By contrast, Neu5, predominantly in BALF, showed higher module scores associated with aging and cell death than Neu1-4 subsets (Fig. 5m and Supplementary Fig. 5f, g). Notably, Neu5 also had lower module scores for 'leukocyte adhesion to vascular endothelial cells', with reduced expression of key genes⁵⁰ (SELL and ITGB2) (Fig. 5n, o). To further define the functional contribution of neutrophil subsets to chronic lung pathology in PASC, we compared inflammatory and fibrotic gene signatures across Neu1–Neu5 populations. Among these, Neu5 exhibited the highest enrichment for gene modules associated with lung fibrosis, TNF-α/NF-κB signaling, and inflammatory responses (Fig. 5p and Supplementary Fig. 5h). Expression of fibrosis-associated genes, including TGFB1, HGF, IL1B, CCL3, CSF1 and CEBPB, was strongly elevated in Neu5 relative to other subsets (Supplementary Fig. 5i). Notably, Neu5 also showed features of cellular senescence and reduced leukocyte adhesion, suggesting impaired clearance and selective retention in BALF. These findings support Neu5 as a functionally distinct, aged neutrophil population that plays a predominant role in driving tissue damage and chronic inflammation in SARS2_PASC. These findings highlight the critical role of neutrophil dynamics in PASC pathogenesis and the therapeutic potential of targeting neutrophil-driven inflammation.

Neutrophil-mediated inflammation has been implicated as a central driver of PASC. Recent studies have identified key proinflammatory mediators, such as S100A8, S100A9 and their heterodimer S100A8/A9 (calprotectin)⁵¹, as well as FPR2, a chemotaxis regulator crucial for neutrophil migration and infiltration at infection sites^52,53. In addition, neutrophil elastase plays a pivotal role in acute lung injury by promoting neutrophil-mediated tissue damage⁵⁴. Based on these findings, we hypothesized that targeting these pathways could attenuate PASC progression.

To evaluate the anti-inflammatory effects of these inhibitors, we measured transcript levels of inflammation-related cytokines (IL1B, NFKB1 and IFNG) in lung and spleen tissues using qRT–PCR. At 15 dpi, 4 days post-treatment, Paquinimod and Sivelestat significantly reduced cytokine expression in the lungs, with Sivelestat showing the greatest reduction, correlating with reduced early mortality (Fig. 6h). This suppression persisted, with cytokine levels remaining low at 30 dpi, 19 days after treatment discontinuation (Fig. 6i). In the spleen, all three drugs reduced IL1B expression at 15 dpi, while Paquinimod uniquely reduced IFNG expression (Supplementary Fig. 6d). By 30 dpi, inflammatory cytokine levels were significantly reduced across all treatment groups in the spleen (Supplementary Fig. 6e). Drug treatment also suppressed neutrophil-associated inflammatory genes (FPR2, S100A9 and ELANE) across tissues, although lung-specific reductions were observed earlier, between 15 dpi and 30 dpi (Supplementary Fig. 6f–i). This difference probably reflects the localized accumulation of neutrophils in the lung, which was absent in the spleen. Among the inhibitors, Sivelestat exhibited the most robust anti-inflammatory effects in both lung and spleen, with suppression of inflammatory gene expression persisting up to 19 days post-treatment. This sustained efficacy underscores Sivelestat's potential as a therapeutic strategy for mitigating PASC, providing a foundation for targeting neutrophil-driven inflammation to alleviate long-term complications of SARS-CoV-2 infection.

Overall, these results highlight that therapeutic strategies targeting neutrophil activity not only reduce inflammation-associated gene expression but also markedly attenuate chronic lung pathology, including fibrosis, induced by SARS-CoV-2 infection. This approach could prove essential for managing lung damage and mitigating long-term complications associated with viral infections.