DMD mRNA and protein expression data were obtained from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/gds) and ArrayExpress (https://www.ebi.ac.uk/arrayexpress/) databases. GSE38417 and E-MTAB-2564 were selected for discovery dataset. A total of 22 samples were collected from GSE38417 (mRNA), including 16 DMD and 6 normal skeletal muscle tissues. A total of 127 DMD and 58 normal plasma samples were collected from E-MTAB-2564 (protein). Then, we downloaded GSE1004 (12 DMD and 12 normal), GSE3307 (10 DMD and 18 normal) and GSE6011 (23 DMD and 14 normal) for validation dataset.

Different intergroup comparisons were performed using the R software (v.3.6.0; https://www.r-project.org/). At first, the mRNA expression matrices of the GSE38417 and protein expression of E-MTAB-2564 datasets were normalized separately. Then, volcano plots and heatmaps of the differentially expressed mRNAs (DEMs) and proteins (DEGs) were plotted by https://www.bioinformatics.com.cn, a free online platform for data analysis and visualization. DEMs with |fold change| ≥ 1.5 and adjusted P-value <0.05 were considered to be statistically differentially expressed. DEPs with |fold change| ≥ 1.2 and P-value <0.05 were considered to be statistically differentially expressed. Moreover, orthogonal partial least squares - discriminant analysis (orthoPLS-DA) was performed to reduce the dimensionality and evaluate the independence of each group through MetaboAnalyst5.0 online tool (16).

To identify modules of highly correlated genes and hub genes, we constructed scale-free co-expression networks for the GSE38417 and E-MTAB-2564 datasets separately through WGCNA (14). First of all, Pearson correlation-based hierarchical clustering was applied to cluster all the genes and samples. Then, soft threshold power value of network construction was selected to make the co-expression network in accordance with a scale-free co-expression network. Third, the adjacency matrices were transformed into a topological overlap matrix for gene modules recognition. In the following step, similar modules were clustered and merged. At last, the module-trait relationship diagrams were performed and the gene list of each module was acquired.

The top two modules, which are most closely related to DMD, were obtained through the module-trait relationship diagrams of the GSE38417 and E-MTAB-2564 datasets. Venn diagram analysis was further used to find overlapping genes between DEGs (DEMs and DEPs) and the above module genes, which were considered as candidate hub genes of DMD and applied for functional enrichment analysis through Metascape software (17).

xCell (18) is a well-established method for calculating the immune and stromal cell composition of “bulk tissue” from microarray or RNA-seq data matrix, which includes 64 stromal and immune cell types. To obtain the differentially expressed cells in DMD, the cell enrich score of stromal and immune cells in the DMD samples were extracted and the difference was performed using Wilcoxon test. Ultimately, Pearson correlation test analysis was implemented to illustrate the relationship between hub genes and differential infiltrated immune and stromal cells, and |r| > 0.7 was set as a strong correlation threshold for further analysis.

SIGNOR 2.0 (https://signor.uniroma2.it/) (19), a public repository that stores manually-annotated causal relationships between proteins and other biologically relevant entities, was performed to explore the signal transduction relationship of the former hub genes. Then, TF-miRNA coregulatory network analysis was implemented on the hub genes through curated regulatory interaction information collected from the RegNetwork repository (20) with Networkanalyst online tools (21).

Intersecting proteins in DEPs and turquoise module from WGCNA were subjected to LASSO regression analyses in 127 DMD and 58 normal samples with glmnet package (22). The ten-fold cross-verification was employed to tune parameter selection. Lambda was used as the minimum partial likelihood deviance. Then, Receiver Operator Characteristic (ROC) curves were plotted to explore the sensitivity and specificity of the intersecting proteins for DMD diagnosis, which was validated in the other four validation datasets, including GSE1004, GSE3307, GSE6011 and GSE38417 datasets. AUC > 0.8 was set as an excellent diagnostic biomarkers threshold for further clinical and animal experiment validation.

To further confirm the results of the proteomic analysis, four hub proteins were validated in 10 DMD patients and 10 healthy controls by ELISAs, including ATP6AP2 (Cat#F15798, YOYOBIO), CTSS (Cat#F08116, YOYOBIO), VIM (Cat#F11800, YOYOBIO), and TNFRSF1B (Cat#F08121, YOYOBIO). Since E-MTAB-2564 used plasma samples for screening biomarkers, we also collected plasma samples of all participants, as well as the characteristics of each patient ( Supplement Table 1 ). Pearson’s correlation analysis was performed to analyze the correlation between the four proteins and the clinical characteristics, including age, white blood cell count (WBC), red blood cell (RBC), platelet count (PLT), C-reactive protein (CRP), alanine amiotransferase (ALT), aspartate aminotransferase (AST), hydroxybutyric dehydrogenase (HBD), creatine kinase (CK), creatine kinase-MB (CK-MB), lactate dehydrogenase (LDH), Motor Function Measure Scale (MFM) score and Mercuri score (23). This study was approved by the Medical Ethics Committee of the first affiliated hospital of Soochow University (No. 2020-145). Informed consent was obtained from all participants.

To further verify the above results, we examined the mRNA expression of the four hub genes of three hindlimb muscles (extensor digitorum longus, EDL; flexor digitorum brevis, FDB; soleus, SOL), comparing 2- and 5-month old WT and mdx mice to investigate degree and time-course of mdx pathology in GSE162455 dataset (55 samples). In addition, mdx mice on C57BL/10ScSnJNju background (mutation in exon 4 of the dystrophin) and wild type mice were bred and used (purchased from GemPharmatech Co. Ltd., n = 6/group). The protocols were approved by the Ethical committee from Soochow University. Experiments were conducted on males at 6–8 weeks of age. Dystrophin protein expression was quantified to evaluate the DMD modeling by immunofluorescence. Proteins were extracted from gastrocnemius muscle. Using SDS-PAGE, the proteins were separated and transferred to PVDF membranes. Incubation with the primary antibodies was carried out at 4°C overnight after blocking with 5% BSA. The secondary antibody was incubated for 60 min. The blots were visualized using the ECL Plus kit (NOBLEBIO, HP5002, Shanghai, China) and exposed to a ChemiDoc MP Imaging System (Bio-Rad, CA, United States). Protein density was calculated using Image Lab software. The variation of protein density was expressed as fold changes compared to the control in the blot. The antibodies used in this study were as follows: CTSS (Cat#A1874, 1:1000), VIM (Cat#A19607, 1:1000), ATP6AP2 (Cat#A6531, 1:1000), GAPDH (Cat#AF7021, 1:1000). Since utrophin is highly related to dystrophin and can substitute for dystrophin’s function (24), we explored the expression levels of these four genes in dystrophin-deficient mice with a transgene expressing high level of full length utrophin.

The “regplot” package was employed to build a nomogram model based on the three hub potential biomarkers. The calibration curve was further used to evaluate the predictive power of the model. Last but not least, decision curve analysis (DCA) and clinical impact curve (CIC) were employed to estimate the clinical value of the model.

The GSE38417 dataset contained the gene expression matrix, including 16 DMD and 6 normal skeletal muscle tissues. The E-MTAB-2564 dataset contained the protein expression profiles from 127 DMD and 58 normal plasma samples. After data preprocessing and differential expression analysis, we identified 4515 DEMs from the GSE38417 dataset (|Fold change| ≥ 1.5 and adjust P value < 0.05); 2652 were up-regulated and 1863 were down-regulated, including MYH8, COL3A1, FZD10, TYRP1, IRX5 and FAM179A ( Figure 2A ). We also identified 56 DEPs from the E-MTAB-2564 dataset (|Fold change| ≥ 1.2 and P value < 0.05); 23 were up-regulated and 32 were down-regulated, including MDH2, MYL3, CA3, RELB, SUCLA2 and VAMP5 ( Figure 2B ). The top 50 DEMs( Figure 2C ) and DEPs ( Figure 2D ) from each dataset were visualized in heatmaps.

To distinguish the significant difference between normal and DMD samples, Orthogonal Partial Least Squares - Discriminant Analysis (orthoPLS-DA) was performed to reduce the dimensionality and evaluate the independence of each group. The results showed that normal samples vs. DMD samples in the two datasets displayed a significant difference ( Figures 3A, C ). In addition, VIP scores were used to distinguish the most important mRNAs and proteins between the groups. The top 10 important features of mRNA were COL3A1, TYRP1, TTC40, COL1A1, CRISPLD1, SPARC, VCAM1, FZD10, COL5A2 and IRX5. The top 10 important features of protein were RELB, SUCLA2, GCC1, CK, POSTN, MDH2, TNFRSF1B, ATP6AP2, H6PD and VAMP5 ( Figures 3B, D ). A total of 25 overlapping genes were detected between the DEMs and DEPs using VIP score ≥ 1( Figure 4A ). GO enrichment analysis showed that these genes were mainly enriched in antigen processing and presentation, collagen−containing extracellular matrix, and endopeptidase inhibitor activity. Pathway enrichment analysis showed that these genes were mainly enriched in adaptive immune system, muscle contraction, cytokine signaling in immune system, and extracellular matrix organization, indicating immune dysfunction in DMD and were used for the subsequent analyses ( Figures 4B, C ).

Genes with similar expression patterns tended to exert similar biological functions. Therefore, we performed WGCNA according to the transcriptome and proteome expression matrix. The hierarchical clustering dendrograms of the samples are shown in Figures 5A, B . Then, we analyzed the soft threshold powers of the network topology and selected β = 5 and 3 as the optimal soft-thresholding parameters ( Figures 5E, F ). We identified 17 modules from the GSE38417 dataset ( Figure 5C ) and 3 modules from the E-MTAB-2564 dataset ( Figure 5D ), and plotted module-trait diagrams to explore the relationships between gene modules and DMD. The module with the highest correlation with DMD was MEturquoise in GSE38417 (r=0.97, p=1e−13) and MEturquoise in E-MTAB-2564 (r=0.54, p= 2e−15) as shown in Figures 5C, D .

To discovery functional genes concurrently involved in DMD across transcriptome and proteome, two methods were applied to achieve the key genes. WGCNA was used to identify DMD-related modules. Differential expression analysis was performed to obtain the most dysregulated genes. In total, we identified 14 functional genes as concurrently involved in transcriptome and proteome, including the turquoise module of the GSE38417 (n = 5441), the turquoise module of the E-MTAB-2564 (n = 247), DEMs of GSE38417 (n = 4515), and DEPs of E-MTAB-2564 (n = 55) ( Figure 6A ). Enrichment analysis showed that these genes were mainly enriched in pathways related to Interleukin-4 and Interleukin-13 signaling, neutrophil degranulation, TNFR2 non-canonical NF-κB pathway and pentose phosphate pathway ( Figure 6B ).

We observed 33 significant differences in the infiltration of cells between DMD and control tissue ( Figure 7A ). In DMD, 24 cell types were up-regulated, including aDC, Adipocytes, Astrocytes, CD4+ memory T−cells, CD8+ Tem, cDC, CMP, DC, Fibroblasts, GMP, HSC, iDC, Keratinocytes, Macrophages, Macrophages M1, Memory B−cells, Mesangial cells, Monocytes, Myocytes, Neurons, Platelets, Preadipocytes, Sebocytes, and Skeletal muscle cell ( Figure 7B ). The nine down-regulated cell types were: CD4+ T−cells, CD8+ naive T−cells, Class−switched memory B−cells, Eosinophils, MSC, Neutrophils, NKT, Pericytes, and Th1 cells ( Figure 7B ).

|r| > 0.7 was set as a strong correlation threshold for further analysis. CTSS was positively correlated with macrophage, DC, macrophage M1 and aDC. VIM was negatively correlated with plasma cells. CD14 was positively correlated with monocytes, DC, and macrophage M1; negatively correlated with D8+ naive T−cells. IFI30 was positively correlated with macrophages. TNFRSF1B was positively correlated with aDC ( Figure 8 ).

NetworkAnalyst online platform is used to generate TF–miRNA co-regulatory network. The interaction between TFs and miRNAs and the selected hub genes was delivered via the TF–miRNA co-regulatory network analysis. This interaction may be the underlying mechanism for the regulation of hub genes expression. The TF–miRNA co-regulatory network is consisted of 70 nodes and 141 edges. A total of 22 miRNAs and 23 TF genes interacted with the validated hub genes, including TFAP2C, MAX, MYC, NFKB1, RELA, hsa-miR-1255a, hsa-miR-130a, hsa-miR-130b, hsa-miR-152, and hsa-miR-17. Figure 9 shows the TF–miRNA co-regulatory network. RELB owned the largest number of neighbors, followed by STAT1, VIM, TNFRSF1B, DES, CD14, CTSS, IFI30, and C3.

As shown in the Figure 10 , we constructed the signaling information network of DEGs based on SIGNOR2.0. STAT1 showed the highest number of interaction, followed by VIM, C3, DES, TNFRSF1B, RELB, CTSS and CD14. Moreover, STAT1 interacted with more complex than other hub genes, including JAK1/STAT1/STAT3, IFNGR2/INFGR1, NfKb-p65/p50, ISGF3 complex. In addition, the analysis also identified the effects of the proliferation, M1 polarization, B cell maturation, and ECM disassembly as potential phenotype, which caused by the alterations in the shortlisted genes. Interestingly, four hub genes (VIM, DES, TNFRSF1B and STAT1) show the interactions through other three genes: AURKB and ROCK1 can directly inhibit VIM and DES; TNFRSF1B can activate STAT1 through increased MAPK14.

To establish a LASSO protein signature for diagnosing DMD, 19 intersecting proteins in DEPs and turquoise module from WGCNA were included for LASSO regression analysis in 127 DMD and 58 normal samples. After 10-fold cross-validation, a LASSO model was established, containing ATP6AP2, CTSS, H6PD, PTCD3, and RELB ( Figure 11A ). To determine which intersecting proteins have the diagnose significance of DMD patients, the ROC analyses were conducted to explore the sensitivity and specificity of the intersecting proteins for DMD diagnosis. The results showed that eight proteins (ATP6AP2, RELB, CTSS, VIM, H6PD, PTCD3, STAT1, and TNFRSF1B) presented the best diagnostic value for differentiating the patients with DMD from healthy controls (AUC > 0.8) ( Figure 11B ). We then validated the intersecting genes in the four validation datasets. And four genes (ATP6AP2, CTSS, VIM, and TNFRSF1B) showed excellent diagnostic performance on discriminating DMD in GSE1004, GSE3307, GSE6011 and GSE38417 datasets (all AUC > 0.8) ( Figure 11C ). This indicated that ATP6AP2, CTSS, VIM, and TNFRSF1B could act as biomarkers to estimate the activity of DMD and verify the effectiveness of the treatment of DMD.

To further increase the reliability of the results, the above proteins (ATP6AP2, CTSS, VIM, and TNFRSF1B) were selected for validation in 10 DMD patients and 10 healthy controls using ELISA. The results showed that ATP6AP2, CTSS, and VIM presented significant higher expression in DMD than healthy controls (P < 0.05) ( Figures 12A–D ). Moreover, the expression of ATP6AP2, CTSS, and VIM in patients with exon 55 mutations was higher than that in patients without exon 55 mutations (P < 0.05) ( Figure 12E ). In addition, LDH was the most related signature between ATP6AP2, CTSS, and VIM expression and characteristics of DMD patients. CRP was the most related signature between TNFRSF1B expression and characteristics of DMD patients ( Supplement Figure 1 ). These observations confirmed that the expression levels of ATP6AP2, CTSS, and VIM are highly specific and can be used as sensitive biomarkers for diagnosing patients with DMD, especially for patients with exon 55 mutations.

To further verify the above results, we first detected the mRNA expression levels of the 4 potential biomarkers in the three mdx mouse muscles, including the extensor digitorum longus (EDL), the flexor digitorum brevis (FDB), and the soleus (SOL). Compared with the sham group, the mRNA expression levels of ATP6AP2, CTSS, VIM, and TNFRSF1B were higher in the three mdx mouse muscles at 2- and 5-month old mice ( Figure 13A ). Since utrophin is highly related to dystrophin and can substitute for dystrophin’s function, we explored the expression levels of these four genes in dystrophin-deficient mice with a transgene expressing high level of full length utrophin. Furthermore, compared with the mdx group, the expression of ATP6AP2, CTSS and VIM recovered in utrophin group ( Figure 13B ). To further verify the results of bioinformatics analysis, we detected the protein expression levels of the 3 potential biomarkers in the muscle tissues. As illustrated in Figures 14A, B , the results showed that the DMD model was successfully established. Compared with the control group, the protein expression levels of CTSS and VIM were significant higher in the mdx group (P < 0.001). The changing trend of ATP6AP2 expression was also consistent with the results of bioinformatics analysis ( Figures 14C–F ). Therefore, ATP6AP2, CTSS and VIM may serve as diagnostic and therapeutic biomarkers.

Then, the nomogram model was constructed to evaluate the probability of DMD based on the 3 significant diagnostic proteins (ATP6AP2, CTSS and VIM) ( Figure 15A ). The calibration curve suggested brilliant agreement among the apparent curve, bias-corrected curve, and ideal curve ( Figure 15B ). DCA and CIC showed that patients could benefit from the nomogram ( Figure 15C, D ). These results indicated that the prediction efficiency of the model was good.