Peripheral blood was obtained from 16 DM-ILD patients and 16 healthy controls. There was no significant difference in the general characteristics between the disease group and the control group (Supplementary Table S1). There was no significant difference in age and clinical data in the DM-ILD patients (Table 1, Table 2). All DM-ILD patients were treated with glucocorticoids and immunosuppressive agents, including prednisone acetate, methylprednisolone, methotrexate, cyclophosphamide, tacrolimus, and cyclosporin A.

Whole blood samples were collected from six DM-ILD patients (ID: LR20L26D X196, LR20L26DX197, LR20L26DX198, LR20L26DX199, LR20L26DX200, and LR20L26DX201) and six healthy controls (ID: LR20L26DX202, LR20L26DX203, LR20L26DX204, LR20L26DX205, LR20L26DX206, and LR20L26DX207) for RNA- seq. Total RNA was isolated and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer’s protocol. The RNA amount and purity of each sample were quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, United States). The RNA integrity was assessed by Bioanalyzer 2100 (Agilent, CA, United States), with a RIN number >7.0, and confirmed by electrophoresis with denaturing agarose gel. Poly(A) RNA is purified from 1 μg total RNA using Dynabeads Oligo (dT)25-61005 (Thermo Fisher, CA, United States) using two rounds of purification. Then the poly(A) RNA was fragmented into small pieces using the Magnesium RNA Fragmentation Module (NEB, Cat. e6150, United States) under 94°C for 5–7 min. Then the cleaved RNA fragments were reverse-transcribed to create the cDNA by SuperScript™ II Reverse Transcriptase (Invitrogen, Cat. 1896649, United States), which were next used to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB, Cat. m0209, United States), RNase H (NEB, cat. m0297, United States), and dUTP Solution (Thermo Fisher, Cat. R0133, United States). The A-base is then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single- or dual-index adapters are ligated to the fragments, and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme (NEB, Cat. m0280, United States) treatment of the U-labeled second-stranded DNA, the ligated products are amplified with PCR by the following conditions: initial denaturation at 95°C for 3 min; eight cycles of denaturation at 98 C for 15 s, annealing at 60°C for 15 s, extension at 72°C for 30 s; and then final extension at 72°C for 5 min. The average insert size for the final cDNA library was 300 ± 50 bp. At last, we performed the 2 × 150 bp paired-end sequencing (PE150) on an Illumina NovaSeq™ 6000 (LC-Bio Technology CO., Ltd. Hangzhou, China) following the vendor’s recommended protocol. The quality of sequencing data was evaluated by FastQC (version 0.11.9). The raw data were filtered by Trimmomatic (version 0.39). Then clean data were mapped to the reference genome of Homo sapiens GRCh38 using HISAT2 (version 2.2.1) with default parameters. Gene counts were calculated by FeatureCount (version 2.0.0). DEGs were identified by R package DESeq2 (version 1.26.0) with | log2 (fold change) | ≥1 and p-value < 0.05. AS events were analyzed by RMATS (version 3.1.1), and differentially expressed AS events were screened with FDR ≤0.05 and |IncLevelDifference| ≥ 0.1 (Shen et al., 2012) (Shen et al., 2014). Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and disease ontology (DO) enrichment analysis were analyzed using R package clusterProfiler. The quantitative analysis of interactive sets between different types of AS was displayed in UpSet plot. The detail of differentially expressed genes (DEGs) and differentially expressed AS events in the chromosome was visualized in the Circos plot.

After DO analysis, DEGs that were enriched into both skin and lung diseases were selected and identified as candidate key genes. Then by overlapping candidate key genes with parent genes of differentially expressed AS, key genes associated with DM-ILD were identified. Key genes and candidate key genes were input into the STRING database to analyze the interactions, and the PPI network was visualized by Cytoscape software (version 3.6.1).

CIBERSORT was applied to assess the infiltration of 22 immune cells in different subgroups, including naive B cells, memory B cells, plasma cells, CD8 T cells, naive CD4 T cells, resting memory CD4 T cells, activated memory CD4 T cells, follicular helper T cells, regulatory T cells (Tregs), gamma delta T cells, resting NK cells, activated NK cells, monocytes, M0 macrophages, M1 macrophages, M2 macrophages, resting dendritic cells, activated dendritic cells, resting mast cells, activated mast cells, eosinophils, and neutrophils. Then the t test was used to compare the difference of immune infiltration between six DM-ILD patients and six healthy controls. The correlations between key genes and differentially distributed immune cells were calculated.

Total RNA of ten DM-ILD patients and ten control subjects were isolated from peripheral blood using TRIzol reagent (15596018, Life Technologies, American). And RNA was reverse-transcribed to cDNA via the SureScript-First Strand cDNA Synthesis Kit (QP057, GeneCopoeia, Guangzhou, China). The expression of PLAUR was detected by Bio-Rad CFX96 (American) using BlazeTaq™ SYBR^® Green qPCR Mix 2.0 (QP033, GeneCopoeia, Guangzhou, China). PCR amplification cycle conditions were 95°C 30 s, 95°C 10 s, 60°C 20 s, and 72°C 30 s for 40 cycles. GAPDH was used as the internal standardized reference. The relative expression level was calculated using the 2^−ΔΔCt method. The specific primers were as follows: PLAUR, 5′-CGAGGTTGT GTGTGGGTTA-3’ (forward) and 5′-GGC​ACT​GTT​CTT​CAG​GGC​T-3’ (reverse); GAPDH, 5′- CGC​TGA​GTA​CGT​CGT​GGA​GTC-3′ (forward) and 5′- GCTGATGAT CTTGAGGCTGTTGTC-3′ (reverse).

All data were analyzed by R (version 4.0.0). Student’s t test was used to compare differences between healthy controls and DM-ILD. Pearson’s correlations between the expressions of key genes and abundance of differentially infiltrated immune cells were calculated. P-values less than or equal to 0.05 were considered statistically significant.

The quality scores, sequence content, and duplication levels of raw sequencing data were evaluated by FastQC and are shown in Supplementary Material S1. To obtain high-quality data for following analysis, low-quality reads were discarded, and the reads contaminated with adaptor sequences were trimmed using Trimmomatic. Thereafter, quality control of clean data was performed by FastQC and is shown in Supplementary Material S2. Statistical results of raw data and clean data are shown in Supplementary Material S3, proving that clean data can be used for following mapping. Then the clean data were mapped to the GRCh38 reference genome, and we found that the percentage of mapped reads was over 79% (Supplementary Material S4), indicating that the clean data meet the needs for following analyses. Next, gene counts were evaluated by FeatureCount (Supplementary Material S5), and the expressions of genes in each sample were calculated. Principal component analysis by R package ggplot2 showed that samples within the same group had good replication (Supplementary Figure S1). A total of 2,018 DEGs were found, including 1,084 upregulated and 934 downregulated genes in DM-ILD patients relative to healthy ones (Figure 1A). The expressions of top 50 upregulated and top 50 downregulated genes in each sample were displayed in the heatmap (Figure 1B).

Thereafter, we investigated the biological function of DEGs in DM-ILD. We found that DEGs were mainly enriched into immune and inflammation biological processes, such as neutrophil activation, neutrophil-mediated immunity, neutrophil degranulation, neutrophil activation involved in immune response, leukocyte migration, humoral immune response, leukocyte chemotaxis and response to the molecule of bacterial origin (Figure 2A), and KEGG pathways, such as cytokine–cytokine receptor interaction, ECM-receptor interaction, IL-17 signaling pathway, NF-kappa B signaling pathway, PI3K-Akt signaling pathway, NOD-like receptor signaling pathway, complement and coagulation cascades, and JAK-STAT signaling pathway (Figure 2B). All of these results indicated that DEGs involved in DM-ILD had a strong relationship with immunity and inflammation. To explore the connection of DEGs and diseases, we further performed DO analysis. We found that DEGs were associated with multiple diseases, such as lung disease, skin disease, arteriosclerotic cardiovascular disease, bacterial infectious disease, obstructive lung disease, and female reproductive system disease (Figure 2C).

By analyzing the transcriptome data, five AS forms, exon skipping (SE), alternative 5 and 3′ splice sites (A5SS and A3SS), mutually exclusive exons (MXE), and intron retention (RI) were identified in DM-ILD patients and healthy ones. After comparison, we obtained 866 differentially expressed AS events from 697 genes, including 97 A3SS in 92 genes, 53 A5SS in 53 genes, 97 MXE in 79 genes, 44 RI in 44 genes, and 575 SE in 494 genes (Figure 3A), indicating that one gene had one or more AS events. Therefore, we used UpSet plot to visualize the intersecting sets of each AS form. As shown in Figure 3B, most genes (91.25%) had only one AS type and 61 genes had two to four types of AS events, including a combination of MXE and SE, A3SS and SE, A5SS and SE, RI and A3SS, MXE and A3SS, RI and A5SS, A5SS and MXE, RI and SE, RI, A5SS, A3SS and SE, MXE, and A3SS and SE. The volcano plots shown in Figure 3C showed the upregulated and downregulated AS events of each AS form. Specifically, there were 54 upregulated and 43 downregulated A3SS, 23 upregulated and 30 downregulated A5SS, 308 upregulated and 267 downregulated SE, 32 upregulated and 12 downregulated RI, and 37 upregulated and 60 downregulated MXE events. Top 50 upregulated and top 50 downregulated AS events in each AS type were displayed in the heatmaps (Figure 3D).

To screen key genes associated with DM-ILD, we extracted 35 common genes from lung disease and skin disease in DO analysis (Figure 4A), namely, BDKRB2, CAMP, COL1A1, F2RL1, F3, MS4A2, CXCR3, ICAM1, IGFBP3, IGHE, IL1A, IL1B, IL1RN, IL6, CXCL8, IL9R, IL10, IL18, CXCL10, KRT18, LAMA5, LAMC2, LEP, MET, CXCL9, MMP1, MMP2, MMP9, PLAUR, PPARG, RNASE2, RNASE3, SOD2, TLR4, and IL1RL1. Using “Dermatomyositis” or “Pulmonary interstitial fibrosis” or “Inflammation” or “Immunization” as key words on the VarElect database, we found that all of these 35 genes were associated with inflammation and immunization. Thus, these 35 genes were identified as candidate key genes. To further explore the role of AS, we overlapped 35 genes with parent genes of differentially expressed AS events and identified PLAUR as the key gene in AS-dependent DM-ILD (Figure 4B). We found that PLAUR had upregulated SE within chr19:43665316–43665459 in DM-ILD (Figure 4C). According to GO and KEGG pathway enrichment analyses, PLAUR genes were found to participate in complement and coagulation cascades, neutrophil-associated immune response, and coagulation-related biological responses, such as neutrophil activation, neutrophil-mediated immunity, neutrophil degranulation, neutrophil activation involved in immune response, blood coagulation, and negative regulation of coagulation. Next, we constructed the PPI network and found that PLAUR had interactions with multiple candidate key genes, including CXCL8, F3, ICAM1, IL1B, IL6, MET, MMP1, MMP2, and MMP9 (Figure 4D). Moreover, we generated a Circos plot to give a full view of DEGs, differentially expressed AS events, and the detail of PLAUR (Figure 5). By RNA-seq, we found that the abundance of PLAUR was higher in DM-ILD samples than the control ones (Figure 6A). Consistent with our sequencing results, RT-qPCR showed that the expression of PLAUR was significantly elevated in the DM-ILD samples (Figure 6B).

DM-ILD was associated with the immune cells. Given the role of immune cells in the immune system, we evaluated the abundance of 22 immune cell types in each sample from DM-ILD patients and healthy ones (Figure 7A). After removing activated memory CD4 T cells, follicular helper T cells, gamma delta T cells, M0 macrophages, M1 macrophages, resting mast cells, and eosinophils with no or few infiltrations in both DM-ILD patients and health ones, the abundance of the remaining 15 immune cell types was compared. We observed that there were significantly more neutrophils and less naive B cells in DM-ILD patients (Figure 7B and Supplementary Figure S2). Given that PLAUR was involved in immune-related biological processes and pathways, we calculated the correlations between PLAUR and neutrophils or naive B cells. Notably, a strong positive correlation (cor = 0.93) was observed between the expression of PLAUR and the abundance of neutrophils (Figure 7C).