A multi-step integrative bioinformatics and experimental validation approach was employed ( Figure 1 ).Transcriptomic data related to CRL and HTS were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), including multiple datasets used for subsequent analysis.The CRL dataset (GSE255848) comprised five subcutaneous adipose tissue samples from CRL patients and five matched healthy controls. The HTS dataset (GSE178411) included 108 samples. For this study, we selected 28 HTS samples and 24 normal controls after excluding non-relevant samples, such as wound healing samples, which were not pertinent to the fibrotic processes under investigation. To validate candidate biomarkers, two additional datasets were analyzed: (1) a scRNA-seq dataset of stromal vascular fraction (SVF) cells from CRL and healthy donors (GSA: HRA000901) (https://ngdc.cncb.ac.cn/gsa-human/browse/HRA000901) (15); and (2) a bulk RNA-seq dataset containing 10 CRL and 10 control skin samples, further stratified by disease duration between short-term (<3 years) and long-term (>14 years) (17).

Differentially expressed genes (DEGs) were independently identified for the CRL (GSE255848) and HTS (GSE178411) datasets using the DESeq2 package in R (v4.4.2). Thresholds were set at adjusted p < 0.05 and |log2FC| > 0.5. DEGs were visualized using volcano plots via the ggvolcano package. Shared DEGs between the two conditions were identified using a Venn diagram tool (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

WGCNA was performed exclusively on the HTS dataset (GSE178411), which includes 52 samples. The CRL dataset (GSE255848) was not suitable for network construction due to its small sample size (n = 10), which is below the generally recommended minimum for robust module detection. The top 25% most variable genes were selected, and a scale-free topology was achieved with soft-threshold power β = 22 and R² = 0.85. Modules with similar expression profiles were merged using a cut height threshold of 0.25, with the minimum module size set to 30 genes. Modules significantly associated with disease phenotype were retained for cross-validation with DEG results.

Shared DEGs(CRL and HTS) and WGCNA(HTS) were mapped to the STRING database (confidence score > 0.4), and the resulting protein–protein interaction (PPI) network was visualized using Cytoscape ((http://www.cytoscape.org). The top five hub genes were selected based on node degree and connectivity.

To support and cross-validate findings from the PPI network analysis, an exploratory machine learning analysis was performed using two algorithms—Random Forest (RF) and Least Absolute Shrinkage and Selection Operator (LASSO). These analyses were conducted on the top five hub genes identified from the PPI network, rather than the full transcriptome, to reduce dimensionality and mitigate overfitting due to limited sample size. RF was implemented using the randomForest package, and LASSO regression via the glmnet package in R. Genes consistently selected by both algorithms were considered supportive fibrosis-associated candidates and are reported in the Supplementary Material for reference.

The xCell algorithm was used to estimate the relative abundance of 64 immune and stromal cell types based on bulk RNA-seq data. Immune profiles of CRL and HTS samples were visualized with ggplot2 to reveal immune microenvironmental alterations.

Functional enrichment analysis was conducted on CRL-associated DEGs, with particular focus on genes expressed in ADSCs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the clusterProfiler package, with adjusted p-values < 0.05 considered statistically significant (18, 19). Enrichment results were displayed via lollipop plots.

All Venn diagrams in this article were generated using an online tool.(https://bioinfogp.cnb.csic.es/tools/venny/index.html).

scRNA-seq data (HRA000901) were processed using Seurat (v4.4.0) and Loupe Browser 8.

Quality control steps included excluding cells with >20% mitochondrial gene content, <500 detected genes, <1,000 UMI counts, or >6,000 genes (to avoid doublets). We applied DoubletFinder (v2.0.3) to identify and remove potential doublets, estimating doublet rates based on the number of loaded cells. Cells with >50% ribosomal gene content were also excluded. To mitigate ambient RNA contamination, we used SoupX (v1.6.2) with default parameters prior to data normalization. After QC, the top 3000 most variable genes were normalized using variance-stabilizing transformation, followed by PCA, clustering (FindNeighbors, FindClusters), and UMAP visualization. Gene expression was explored using FeaturePlot. Intercellular communication networks were inferred with CellChat (v1.5.0), focusing on ligand–receptor interactions relevant to fibrosis and inflammation. AUCell (v1.22.0) was used to assign fibrosis and inflammatory scores to each single cell based on curated gene sets.

A total of 11 female patients diagnosed with CRL were recruited from the Department of Burns and Plastic Surgery, Affiliated Hospital of Zunyi Medical University between January 1, 2024, and January 1, 2025, for validation. Inclusion criteria were: age between 35 and 70 years, confirmed stage II–III lower limb CRL, disease duration of at least 2 years, and no concurrent systemic fibrotic or autoimmune conditions. All patients were in the chronic phase of CRL-associated fibrosis. This study was approved by the Institutional Ethics Committee of the Affiliated Hospital of Zunyi Medical University, and written informed consent was obtained from all participants. Tissue samples were collected as discarded material during clinically indicated surgical procedures. No additional intervention beyond routine care was performed. Ethical principles and patient autonomy were fully respected during recruitment.

Paraffin-embedded tissue sections (4 μm thick) were stained with Masson’s trichrome to evaluate collagen deposition and fibrotic remodeling. Sections were imaged using bright-field microscopy. For each sample, three representative high-power fields were selected, and the percentage of collagen fibers relative to the total tissue area was calculated and averaged.

Human skin tissue samples were collected from consenting individuals. Total proteins were extracted using RIPA lysis buffer and quantified via the BCA assay. Equal amounts of protein were subjected to SDS-PAGE, followed by transfer onto PVDF membranes. Membranes were incubated with primary antibodies against ASPN, TGF-β1, and GAPDH, then with HRP-conjugated secondary antibodies. Signals were visualized using enhanced chemiluminescence (ECL) and analyzed with ImageJ software.

Statistical analyses were conducted using R software (version 4.4.2). Differences between groups were evaluated using unpaired two-tailed Student’s t-tests. Correlation analyses were performed using either Pearson or Spearman correlation coefficients, depending on data distribution. A p-value less than 0.05 was regarded as statistically significant.

To investigate fibrosis-associated transcriptional changes in CRL and HTS, differential expression analysis was first performed. In the GSE255848 dataset (CRL), 963 DEGs were identified, comprising 613 downregulated and 350 upregulated genes ( Figures 2A, B ). Similarly, 9952 DEGs were obtained from the GSE178411 dataset (HTS), including 4390 downregulated and 5562 upregulated genes.

To further uncover disease-relevant gene modules in HTS, WGCNA was performed using the GSE178411 dataset. A soft-thresholding power of β = 22 was chosen to achieve a scale-free network topology ( Figure 2C ), resulting in the identification of 15 distinct gene co-expression modules ( Figure 2D ). Among these, the blue module exhibited the strongest association with HTS status (|r| = 0.84, p = 3e−15) and encompassed 2,272 genes ( Figures 2E, F ), which were selected for subsequent analyses.

We next identified overlapping fibrosis-related genes by comparing three sets: CRL DEGs (963 genes), HTS DEGs (9952 genes), and genes in the blue module from WGCNA (2272 genes). A total of 154 shared genes were identified using an online Venn diagram tool ( Figure 2G ), representing potential mediators of common fibrotic pathways in both conditions.

The 154 shared genes were input into the STRING database to generate a PPI network, which was visualized in Cytoscape ( Figure 3A ). Using the MCODE plugin, five highly interconnected hub genes were identified: POSTN, ASPN, FMOD, COL11A1, and MFAP5 ( Figure 3B ). Further expression analysis revealed that only four of these genes(POSTN, ASPN, FMOD, and MFAP5) exhibited statistically significant differential expression between lymphedema and control tissues ( Figure 3C ), indicating their potential biological relevance in the fibrosis pathogenesis of lymphedema.

To further support hub gene prioritization, we applied RF and LASSO models to the five genes identified from the PPI network. LASSO identified ASPN, FMOD, and MFAP5 as candidate predictors in GSE255848. RF feature importance ranking confirmed these same three genes. Venn analysis showed complete overlap between both models. Although based on a limited dataset, this exploratory analysis yielded results consistent with the core PPI findings. Notably, ASPN showed the most robust upregulation in chronic lymphedema tissue. Functionally, ASPN modulates ECM remodeling and TGF-β signaling; FMOD is involved in collagen fibrillogenesis; and MFAP5 contributes to microfibril integrity. These findings are included in Supplementary Figure S1 , and should be interpreted as supportive observations rather than conclusive evidence, due to the small sample size.

Immune microenvironment profiling via the xCell algorithm revealed distinct immune cell infiltration patterns between disease and control groups in both datasets ( Figures 4A, B ). Correlation analysis showed that POSTN, ASPN, FMOD, and MFAP5 were most strongly associated with preadipocyte infiltration in both CRL and HTS samples ( Figures 4C, D ), suggesting a link between mesenchymal transition and immune remodeling in fibrosis.

To assess the biomarker potential of candidate genes, we analyzed an independent bulk RNA-seq dataset from CRL patients stratified by disease duration (<3 years vs. >14 years). Among all candidates(POSTN, ASPN, FMOD, MFAP5), ASPN was the only gene consistently overlapping with differentially expressed genes (DEGs) in long-term cases across multiple comparisons ( Figure 5A ), indicating a strong association with advanced fibrotic stages. To elucidate the cellular context of ASPN expression, we conducted a focused analysis on scRNA-seq data from CRL tissue. ASPN expression was enriched in a specific subset of stromal cells ( Figures 5B, C ). Notably, ASPN expression was also observed in a subset of pericytes, suggesting that these cells may also participate in the fibrotic microenvironment in CRL. Violin plots further demonstrated that ASPN expression was largely confined to ADSCs, with negligible expression in immune and epithelial compartments ( Figure 5D ). GO and keg enrichment of ADSCs cells revealed associations with key fibrotic pathways, including Wnt signaling, TGF-β response, collagen fibril organization, and mesenchymal cell proliferation ( Figures 5E, F ). Additionally, AUCell hallmark gene set analysis showed that ADSCs were enriched for EMT, TGF-β, IL6–JAK–STAT3, and hypoxia signaling pathways ( Figure 5G ), underscoring their active role in extracellular matrix remodeling and inflammation-associated fibrosis. Together, these findings establish ASPN as a ADSCs-enriched marker of a distinct, pro-fibrotic stromal subpopulation, supporting its utility as both a biomarker and a potential therapeutic target in late-stage CRL fibrosis.

To elucidate the role of ASPN in CRL fibrosis, we performed a combined histological and molecular analysis of skin tissue samples from 11 CRL patients with varying disease durations and degrees of fibrosis ( Supplementary Material 1 ). Western blot analysis revealed a increase in ASPN protein expression correlated with longer disease duration, more severe histological fibrosis, and elevated TGF-β1 levels ( Figure 6 , Supplementary Material 2 ). Shapiro–Wilk tests confirmed normal distribution for all variables (Masson staining, disease duration, ASPN, and TGF-β1; all p > 0.05), supporting the use of Pearson correlation. Quantitative image analysis showed a positive correlation between collagen content and disease duration (r = 0.76, p = 0.0032), indicating progressive fibrotic remodeling over time. Importantly, expression levels of ASPN and TGF-β1 were positively correlated not only with disease duration (r = 0.89, p = 0.0001 and r = 0.55, p = 0.033, respectively) but also with the extent of collagen deposition observed in Masson staining (r = 0.76, p = 0.0035 and r = 0.77, p = 0.0026, respectively). In summary, these findings indicate that ASPN expression is associated with TGF-β1 levels and may have potential as a biomarker for disease progression and the severity of fibrosis.

CellChat analysis revealed enhanced intercellular communication in CRL tissues relative to controls ( Figure 7A ). Notably, TNF and IFN signaling pathways were upregulated and primarily derived from immune cell, with ADSCs acting as the main recipient cells ( Figures 7B, C ), suggesting that chronic inflammation may drive ADSC activation and downstream ECM remodeling. In the LSVF(lymphedema stromal vascular fraction) group, this communication network was markedly amplified. Especially, T cells exhibited increased outgoing signaling, while ADSCs received a significantly greater proportion of incoming signals compared to the NSVF(normal stromal vascular fraction) group. Several pro-inflammatory and pro-fibrotic pathways—including CXCL, OSM, IFN, and TNF—were notably enriched in LSVF, reflecting a shift toward sustained fibrotic activation. The elevated signaling from multiple immune subsets to ADSCs likely contributes to persistent ECM deposition and progressive tissue remodeling, characterizing advanced fibrotic stages.