Single-cell RNA sequencing analysis was conducted on PBMC samples from three ASS-ILD treatment-naive patients and three healthy individuals recruited at Shantou Central Hospital. Healthy individuals were selected to be sex-, ethnicity-, and age-matched. Informed consent was obtained from all the subjects. The diagnosis of ASS-ILD was made by a multidisciplinary team comprising an expert rheumatologist and two experienced radiologists specializing in chest CT. ASS was diagnosed in accordance with the criteria proposed by Solomon et al. (8). Patients with other identifiable causes of ILD, including those with medication-related lung injury, malignancy, or environmental and occupational exposures, were excluded from the study. Supplementary Table 1 provides a detailed description of the clinical presentation and laboratory characteristics of the patients, while Supplementary Figure 1 offers a visual representation of the CT images of the patients. All samples were collected in accordance with the ethical requirements and regulations of the Ethics Committee of Shantou Central Hospital. Informed consent was obtained from all the subjects, and the studies were conducted under approval (approval number:〔2022〕KY-006).

Furthermore, a single-cell RNA sequencing dataset [GSE190510 (9)] from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) was incorporated, comprising eight PBMC samples from five ASS-ILD patients and three healthy individuals.

Single-cell sequencing was conducted via the 10x Genomics platform (1). Preparation of single-cell suspensions: The isolation of peripheral blood mononuclear cells (PBMCs) from undiluted human blood was conducted via Histopaque solution (Sigma–Aldrich, St. Louis, MO). A total of 10 µl of the suspension was counted under an inverted microscope with a hemocytometer. The number of live cells was determined via the Trypan blue method (2). Construction of single-cell libraries: Chromium Single-cell 3’ Reagent v3 kits were used to prepare barcoded single-cell RNA sequencing libraries in accordance with the manufacturer’s instructions. The isolated PBMCs were encapsulated via microfluidics technology and barcoded with a unique molecular identifier. cDNA was prepared in accordance with the manufacturer’s specifications (3). Single-cell RNA sequencing library preparation and sequencing: cDNA libraries were sequenced on an Illumina HiSeq PE150 system (4). Raw data processing and quality control: The data were demultiplexed via Cell Ranger software (version 3.1.0), which generated FASTQ files, which were aligned to a human reference genome (GrCh38). The Cell Ranger software generated a unique raw molecular identifier count matrix, which was subsequently converted into a Seurat object via the R package Seurat (version 5.0.1). Cells with doublets and low quality were removed based on the number of unique molecular identifiers (UMIs) and the proportion of mitochondrial gene expression. The genes were filtered on the basis of the number of cells in which they were expressed, and the cells were filtered based on the number of genes expressed in them. The data were subsequently normalized by log normalization, and the top 2000 highly variable genes were selected based on mean expression and variance. All genes were scaled via the ScaleData function, and principal component analysis (PCA) downscaling was performed. The cells were subsequently clustered via the FindNeighbors and FindClusters functions to obtain cell subgroups, and the cells were subsequently annotated. Batch correction was performed via the Harmony algorithm.

The identification of differentially expressed genes and specific marker genes for each cellular subset was achieved by employing Seurat’s FindAllMarkers() function with the parameter “test.use = wilcox” by default under the RNA assay. The definition of each cell subcluster was based on the expression of canonical markers.

To identify the genes that were upregulated in different cell types or different disease states, the FindMarker function (Logfc.threshold = 0.25, Wilcoxon test) in Seurat was employed. To investigate the biological functions and pathways associated with the differentially expressed genes (DEGs, log2-fold change (FC) > 0.25, adjusted p value < 0.05), Gene Ontology Biological Process (GOBP) (10) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) (11) pathway analyses were conducted via the clusterProfiler package (version 4.8.3) (12).

The DEGs of each cluster were then used as input to generate ISGs via the Interferome database (https://interferome.org/interferome/home.jspx). This process yielded 235 ISG gene sets ( Supplementary Table 2 ), which were subsequently used for ISG scoring via the AUCell R package (version 1.24.0) (13). The ISG set was then employed as the input data for the calculation of the area under the curve (AUC) value. The AUC values were then utilized to construct gene expression rankings for each cell. The AUC provides an estimate of the proportion of genes within the gene set that are highly expressed in each cell. The number of expressed genes in a cell was positively correlated with the AUC value. Consequently, cells that express a greater number of genes from the gene set will have higher AUC values than cells that express fewer genes. The function “AUCell_exploreThresholds” was employed to ascertain the threshold that could be utilized to consider the present gene set active. The cell clustering UMAP embedding was subsequently colored according to the AUC score of each cell, thereby indicating which cell clusters were active in the ISG gene set.

Human gene sets from GOBP were retrieved via the msigdbr package (version 7.5.1). We subsequently applied GSVA (14) with the GSVA package (version 1.50.0) to assign pathway activity estimates to individual cells.

The cell–cell interactions between different cell types were evaluated via CellChat (version 1.6.0, R package). CellChat utilizes gene expression data as the fundamental input to model the probability of cell-to-cell communication by integrating gene expression data with an existing database comprising known interactions between signaling ligands, receptors, and their cofactors (15). Normalized count data from each condition were used to create a CellChat object, and the recommended preprocessing functions for the analysis of individual datasets were applied with default parameters.

Gene set enrichment analysis (GSEA) of monocytes from ASS-ILD patients and HCs was performed via the clusterProfiler (version 4.12.6) R software package to analyze the potential biological pathways of monocytes in ASS-ILD. Permutations were set to 10,000 to obtain normalized enrichment scores (NESs) in GSEA. Gene sets with an adjusted P value <0.05 were considered to be significantly enriched. The enrichplot (version 1.24.2) and ggplot2 (version 3.5.1) R packages were employed to display the enrichment results.

Transcription factor (TF) activity was inferred for monocytes via a univariate linear model (ULM) in the R package decoupleR (version 2.12.0) (16), with CollecTRI (Collection of Transcriptional Regulatory Interactions) serving as the reference. CollecTRI is a comprehensive resource that contains a curated collection of TFs and their transcriptional targets, compiled from 12 different resources (17).

We additionally collected PBMC samples from 5 ASS-ILD patients and 5 matched healthy controls. RNA was extracted using the Trizol method, and cDNA was synthesized using Takara’s reverse transcription kit (cat#RR092A) from Takara Bio Inc. (Shiga, Japan). Semi-quantitative PCR was performed using MCE’s SYBR Green qPCR Master Mix (cat# HY-K0501) from MCE (Shanghai, China), and the results were analyzed using the 2−ΔΔCt method with β-actin as the reference gene for normalization. The primers for each gene are as follows:

Additionally, we incorporated an external bulk RNA-seq dataset (GSE220915 (18)), which includes muscle biopsy samples from 18 ASS patients and 33 healthy controls, to validate the expression of key genes.

All statistical analyses were performed using R software (version 4.3.2) or GraphPad Prism (version 10.4.0). Normally distributed data are expressed as the mean ± standard deviation and compared using unpaired t-tests. Non-normally distributed data are presented as median (interquartile range) and analyzed using the Mann-Whitney U test. A p value of less than 0.05 was considered statistically significant.

To elucidate the immunopathological landscape of ASS-ILD, we first analyzed a discovery scRNA-seq cohort comprising peripheral blood mononuclear cells (PBMCs) from three treatment-naive ASS-ILD patients and three healthy controls (HCs), yielding a total of 67,421 cells. Unsupervised clustering followed by UMAP visualization ( Figures 1A, B ) identified eight major immune cell types: CD4⁺ T cells, CD8⁺ T cells, monocytes, natural killer (NK) cells, B cells, low-density granulocytes (LDGs), megakaryocytes (MKs), and dendritic cells (DCs), confirmed by canonical marker gene expression ( Figure 1C ).

To assess interferon activity, AUCell scoring based on 235 curated interferon-stimulated genes (ISGs) from the Interferome database revealed significantly elevated ISG scores in monocytes compared with other immune subsets ( Figures 1D, E ). Gene ontology (GO) and KEGG pathway enrichment further demonstrated that monocytes were enriched for pathways associated with innate immune activation (e.g., IL-6 and TNF production, myeloid cell activation), lysosomal processing, mycobacterial infection, and autoimmune inflammation ( Figures 1F, G ). These findings indicate that monocytes are a key site of interferon-driven immune activation in ASS-ILD.

To enhance statistical power, we integrated our dataset with an external cohort (GSE190510), yielding a total of 126,026 cells for analysis. Consistent with our discovery findings, monocytes displayed the highest ISG activity across all immune lineages in the integrated dataset ( Figures 2A–C ). Subsequent subclustering identified six transcriptionally distinct monocyte subsets (Mono0–Mono5) ( Figure 2D ).

Classical monocytes (Mono0/1) expressed CD14 along with genes related to leukocyte migration (e.g., GRN, LGALS2) and apoptotic recognition (e.g., CD36, TLR2). Mono2 was characterized by upregulation of IFNG, STAT4, and RORA, suggesting a proinflammatory trajectory involving Th1/Th17 polarization and NK cell activation. Mono3 showed proliferative signatures via ribosomal and translational machinery (e.g., RPL5, EEF1B2). Nonclassical CD16⁺ monocytes (Mono4/5) displayed IFN-polarized features: Mono4 was enriched in antigen presentation (e.g., ISG15, HLA), while Mono5 showed IFN-adaptive potential through IFNAR1/2 overexpression. Functional heatmap analysis revealed distinct gene expression and pathway enrichment across subsets, suggesting specialized roles in ASS-ILD pathogenesis ( Figures 2E, F ).

To delineate the developmental trajectory of monocyte subsets, we applied pseudotime inference using the Monocle algorithm. The trajectory branched from Mono0 toward two terminal fates—Mono1 and Mono2. Cells from ASS-ILD patients showed a preferential trajectory toward Mono2, while those from HCs were more diffusely distributed toward Mono0 and Mono1 ( Figure 3A, C ). Notably, ISG scores increased along the pseudotime axis ( Figure 3D ), supporting progressive interferon polarization.

Monocyte subset density across pseudotime revealed early-stage enrichment of Mono0, Mono3, and Mono4, with later stages dominated by Mono1, Mono2, and Mono5 ( Figure 3E ). Differential gene expression over pseudotime (adjusted p < 0.0001, logFC > 1) clustered into three groups enriched for leukocyte proliferation/migration, innate immunity, and adaptive immunity, respectively ( Figure 3F ). These results suggest that IFNG⁺ monocytes (Mono2) accumulate at later pseudotime stages and may contribute to the transition from basal to inflammatory immune phenotypes in ASS-ILD.

To explore intercellular signaling dynamics, we performed CellChat analysis to compare ligand–receptor interactions between ASS-ILD patients and HCs. Although global communication strength was reduced in ASS-ILD ( Figure 4A ), there was selective upregulation of type II interferon (IFN-γ) and tumor necrosis factor (TNF) signaling pathways ( Figure 4B ). Directional analysis revealed that Mono2, NK cells, and CD8⁺ T cells were dominant sources of IFN-γ and TNF signals, targeting monocyte subsets (Mono1/2/4/5) ( Figures 4C, D ).

Network visualization showed that IFN-γ signaling from Mono2, NK, and CD8⁺ T cells converged on downstream monocyte subsets, while TNF signals from Mono1, Mono2, and Mono5 regulated NK cells, CD8⁺ T cells, DCs, and other monocytes ( Figures 4E, F ). Ligand–receptor profiling indicated comparable receptor levels between groups but markedly elevated TNF-α and IFN-γ ligand expression in ASS-ILD ( Figure 4G ). These data reveal that IFN-γ and TNF signaling networks, driven in part by the Mono2 subset, amplify monocyte–immune crosstalk and may contribute to sustained inflammatory circuits in ASS-ILD.

To identify regulatory drivers of monocyte subset behavior, we employed decoupleR-based inference using CollecTRI TF-target interaction data. Six transcription factors—ETV5, IRF5, IRF7, RORB, RORC, and SMAD1—showed the highest inferred activity in ASS-ILD monocytes ( Figure 5A ). These TFs are linked to immune polarization, interferon signaling, and fibrosis: ETV5 promotes Th1/Th17 responses via STAT3/STAT4 activation and IFN-γ/IL-17A production (19); IRF5 and IRF7 are canonical regulators of IFN pathways and macrophage activation (20, 21); RORB and RORC regulate autophagy and Th17 cell differentiation (22–24); SMAD1 mediates TGF-β signaling and contributes to fibrotic remodeling (25).

Gene set enrichment analysis (GSEA) of these TF targets showed enrichment in IL-17, T cell receptor, and inflammatory bowel disease pathways ( Figure 5B ). Notably, Mono2 exhibited the highest TF activity among monocyte subsets ( Figure 5C ). Pathway activity profiling revealed that Mono1, Mono2, and Mono5 in ASS-ILD were enriched in innate and adaptive immunity-related pathways, including IL-17, JAK-STAT, NF-κB, TGF-β, and TNF signaling ( Figure 5D ). These results position Mono2 as a transcriptionally primed subset with proinflammatory and profibrotic potential.

To further validate our findings, we assessed the expression of key genes, including IFNG, TNF, IRF7, and NFKBIZ. The qRT-PCR results revealed that the expression of these genes was significantly higher in the PBMCs of ASS-ILD patients compared to healthy controls ( Figure 6A ). Additionally, we incorporated a bulk RNA-seq dataset from muscle biopsy samples, which included 13 ASS patients and 33 healthy controls. Similar to the findings in PBMCs, the expression levels of these four genes were elevated in the muscle tissue of patients, compared to healthy controls ( Figure 6B ). This suggests that these genes are highly expressed not only in peripheral blood but also in affected tissues, further confirming the reliability of our results.