This study was approved by the Ethics Committee of Shanghai Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine (Approval ID: SH9H-2021-T100-2). All participants provided written informed consent in accordance with the Declaration of Helsinki (version 2002). From January 2017 to December 2018, oral mucosal tissue samples were collected from patients diagnosed with oral lichen planus (OLP) or oral lichenoid lesions (OLL) and from healthy volunteers undergoing orthognathic surgery. Diagnoses were established based on clinical and histopathological criteria as previously described (1, 31–34). OLL cases were strictly defined as lesions with atypical clinicopathologic features after excluding known triggers, including drug-induced, contact-induced, and GVHD-associated lesions (35). Inclusion and exclusion criteria followed our previously reported standards (30). The workflow is presented in Figure 1 .

Formalin-fixed, paraffin-embedded (FFPE) oral mucosal tissue sections from three patients (two with non-erosive and one with erosive OLP/OLL) and three healthy controls were analyzed using the 10x Genomics Visium HD platform. Tissue sections (10 μm) were deparaffinized, permeabilized, and stained with hematoxylin and eosin (H&E) for spatial alignment, following the Visium HD FFPE Tissue Preparation Handbook (CG000684, 10x Genomics). The CytAssist system was used to position the tissue onto Visium HD slides, which contain a high-density array of 2 × 2 μm spatially barcoded capture spots. RNA quality of FFPE tissue was assessed using the DV200 metric, with a minimum of 30% of RNA fragments >200 nucleotides required for high-quality RNA. Probe hybridization, ligation, extension, and library construction were performed according to the Visium HD Spatial Gene Expression Reagent Kits User Guide (CG000685), and libraries were sequenced on an Illumina NovaSeq X Plus-25B platform with paired-end reads.

The raw sequencing data were processed using the Space Ranger pipeline (version 4.0.1). Initially, reads were trimmed to remove template switch oligo (TSO) sequences from the 5’ end and polyA tails from the 3’ end. Barcode correction was performed by matching observed barcodes to a predefined inclusion list, allowing up to four permissible edits. The reads were then aligned to the GRCh38 human reference genome (GENCODE v38 annotation) using the STAR aligner (version 2.7.1b), with high sensitivity to spliced alignments. Alignment quality was assessed by classifying reads as exonic, intronic, or intergenic based on their overlap with annotated genomic features. A filtered feature-barcode matrix was generated by counting unique molecular identifiers (UMIs) for each gene, ensuring the inclusion of only high-quality reads in downstream analysis. Custom quality control filters were applied to exclude low-quality spots, specifically those with fewer than 10 detected features (nFeature_Spatial) or total counts (nCount_Spatial) below 10. These thresholds were chosen to eliminate spots with insufficient gene expression data, thereby reducing noise and improving result accuracy. The data were normalized and variance-stabilized using SCTransform to account for both technical noise and biological variability. Dimensionality reduction was performed using UMAP (Uniform Manifold Approximation and Projection) on the highly variable genes to capture the most informative features of the dataset. Unsupervised clustering was then applied to identify spatially distinct transcriptional domains. Finally, marker gene expression was mapped to reconstruct the local immune architecture within the tissue. For downstream analysis, the data were initially generated at a 2 μm resolution, with additional binning resolutions of 8 μm and 16 μm. To balance cell resolution and the mean transcripts per bin, the 8 × 8 μm bin resolution was selected for further analysis. This resolution offers an optimal compromise between spatial resolution and noise reduction, providing sufficient detail while maintaining data quality.

Spatial transcriptomics gene signature scoring: Tertiary lymphoid structure (TLS) gene lists were compiled from previous studies (24). Gene signature scores for each Visium slide were calculated using the AddModuleScore function in Seurat (v4.4.1) with default settings. Briefly, AddModuleScore computes the average expression of the genes in the set for each spatial spot and subtracts the aggregated expression of control feature sets, which are randomly sampled from bins of genes with similar average expression to control for expression-level biases. The resulting scores represent relative TLS-associated transcriptional activity per spot and were visualized using SpatialFeaturePlot in Seurat, where the color scale corresponds to these relative values.

Oral mucosal biopsies were collected under local anesthesia and bisected for histopathological evaluation. Total RNA was extracted using TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s protocol. RNA purity and concentration were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Scientific, USA), and RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

RNA-seq libraries were constructed using the VAHTS Universal V10 RNA-seq Library Prep Kit (Premixed Version) following the manufacturer’s instructions. Libraries were then sequenced on the Illumina NovaSeq 6000 platform, generating 150 bp paired-end reads. Raw FASTQ reads were processed with fastp to remove low-quality sequences, yielding clean reads for subsequent analyses.

The clean reads were aligned to the reference genome using HISAT2. Gene-level read counts were obtained using HTSeq-count. Bulk RNA-seq was performed on samples from 16 OLP patients and 10 healthy controls. To expand the dataset, two publicly available transcriptomic datasets (GSE213346 and GSE204663) were retrieved from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) and integrated with the in-house cohort.

Raw counts were normalized to Transcripts Per Million (TPM), and batch effects were corrected using the ComBat function in the sva R package. Principal component analysis (PCA) was performed using R (version 4.4.1) to assess biological variation and sample clustering.

A curated set of 39 tertiary lymphoid structure (TLS)-associated genes was compiled based on published literature (24). These genes include chemokines (CCL2/3/4/5/8/18/19/21, CXCL9/10/11/13), T follicular helper (TFH) cell markers (CXCL13, CD200, FBLN7, ICOS, SGPP2, SH2D1A, TIGIT, PDCD1), genes shared by T helper 1 (TH1) and B cell lineages (CD4, CCR5, CXCR3, CSF2, IGSF6, IL2RA, CD38, CD40, CD5, MS4A1, SDC1, GFI1, IL1R1, IL1R2, IL10, CCL20, IRF4, TRAF6, STAT5A), and TNFRSF17 as a canonical plasma cell marker. To characterize immune heterogeneity in OLP, consensus clustering was performed using the intersection of differentially expressed genes (DEGs) and TLS signature genes (TSGs). Unsupervised clustering was conducted using the R package ConsensusClusterPlus, with the number of clusters (k) explored from 2 to 7. To ensure clustering stability, 1,000 bootstrap resampling iterations were performed. In each iteration, 80% of features (Pltem = 0.8) were randomly selected, while all genes were eligible for inclusion (pFeature = 1.0). Hierarchical clustering (clusterAlg = “hc”) with Pearson correlation distance (distance = “pearson”) was applied. The optimal number of clusters was determined based on cumulative distribution function (CDF) plots and delta area analysis, aiming for maximal intra-cluster homogeneity and inter-cluster separation. Clustering consistency was visualized using the ggplot2 package.

Differentially expressed genes (DEGs) were identified using the DESeq2 package in R. Gene symbols were converted to Entrez IDs using the org.Hs.eg.db package. Genes with an adjusted P-value (Padj) < 0.05 and |log²FoldChange| > 1 were considered statistically significant and visualized using a volcano plot to highlight robust transcriptional changes.

To explore the biological functions and pathways associated with the observed gene expression changes, Gene Ontology (GO) enrichment analysis was performed using the ClueGO plugin in Cytoscape (36). For enrichment input, a slightly broader threshold of nominal P-value < 0.05 and |log²FoldChange| > 1 was applied to include additional genes and minimize false negatives. Functional enrichment was assessed using a right-sided hypergeometric test, with P-values corrected for multiple testing using the Benjamini-Hochberg method (36). GO terms were grouped based on kappa statistics (threshold = 0.4) into functionally related categories. Enriched pathways with adjusted P-value < 0.05 were considered statistically significant.

In summary, more stringent criteria were applied for DEG visualization, whereas a more inclusive threshold was used for enrichment analysis to balance specificity and sensitivity across different analytical purposes.

An independent retrospective clinical cohort of 26 OLP and 48 OLL patients was enrolled to validate spatial immune cell infiltration. Disease severity was assessed using the Reticular-Erythema-Ulcerative (REU) scoring system (37), with consistent diagnostic and scoring protocols applied across cases to reduce inter-observer bias.

The formalin-fixed paraffin-embedded (FFPE) tissue sections were deparaffinized in xylene and rehydrated through a graded ethanol series. Endogenous peroxidase activity was quenched by incubating the sections with 3% hydrogen peroxide in methanol for 5 minutes. Antigen retrieval was performed by heating the sections in tris-EDTA buffer (pH 8.0) in a 100°C water bath for 30 minutes, followed by cooling at room temperature and rinsing with tap water. The sections were then incubated with a blocking solution (P0102, Beyotime, Hangzhou, China) at room temperature for 1 hour to reduce nonspecific binding. Following blocking, the slides were incubated overnight at 4°C with primary antibodies: anti-CD20 (1:50, ab78237, Abcam, Cambridge, UK) and anti-CD38 (1:500, ab108403, Abcam, Cambridge, UK). After washing with phosphate-buffered saline (PBS), the sections were incubated with the secondary antibody, Goat Anti-Rabbit IgG H&L (HRP) (1:2000, ab205718, Abcam, Cambridge, UK), for 1 hour at room temperature. Detection was carried out using diaminobenzidine (DAB) (Celnovte, MD, USA) as a chromogen, followed by counterstaining with hematoxylin for 1 minute. The slides were then dehydrated through a graded ethanol series, cleared in xylene, and mounted with a coverslip for microscopic examination.

For semiquantitative evaluation, five representative fields per slide within the inflammatory zone were selected at 400× magnification. Immunostaining was graded as follows: Grade 3 (> 2/3 positive area or follicle-like aggregates), Grade 2 (1/3 to 2/3 positive area), Grade 1 (< 1/3 or focal positivity), and Grade 0 (no detectable staining) (30).

In addition, TLSs were assessed in the same sections using the criteria established in our previous study (29), where positivity was defined by the presence of high endothelial venules (HEVs, MECA-79) and CD3+ T-cell clusters surrounding CD20+ B-cell aggregates. All sections were independently reviewed by two experienced pathologists. In case of discrepancy, consensus was reached via joint discussion.

Continuous variables are presented as mean ± standard deviation (SD), while categorical variables are expressed as counts and percentages. For continuous variables, Welch’s t-test was applied for two-group comparisons, and one-way ANOVA followed by Tukey’s post-hoc test was used for comparisons among three or more groups. Categorical variables were compared using Fisher’s exact test or Chi-squared test, as appropriate. All statistical analyses were conducted in R software (version 4.2.2), with two-sided P < 0.05 considered statistically significant.

A curated panel of 39 tertiary lymphoid structure (TLS) signature genes (TSGs), encompassing chemokines, T follicular helper cell markers, and genes shared by Th1 and B cell lineages, was selected based on previous studies (24). Spatial transcriptomic analysis of oral mucosal tissues (n = 6; 3 OLP/OLL patients and 3 healthy controls) revealed increased expression of TLS-related genes in OLP/OLL lesions, with a stronger signal observed in erosive cases ( Figure 2A ). TLS-like aggregates were visualized based on the spatial co-localization of canonical TLS markers, including CD3E, FDCSP, AICDA, BCL6, MS4A1, and CD38 ( Figure 2B ), supporting the presence of organized immune niches within inflamed tissues.

To further characterize immune heterogeneity in OLP, we integrated three transcriptomic datasets comprising 16 in-house OLP samples, 10 healthy controls, and two public datasets (GSE204663: 19 OLP; GSE213346: 40 OLP). Batch effects were corrected using the ComBat algorithm, and successful integration was confirmed by principal component analysis (PCA) ( Figure 2C ). Differential expression analysis comparing OLP (n = 75) and healthy controls (n = 10) identified 2,350 significantly dysregulated genes ( Figure 2D ). By intersecting these differentially expressed genes (DEGs) with the 39 TSGs, we identified 26 transcripts ( Figure 2E ) with consistent upregulation in the OLP group ( Figure 2F ).

To define TLS-related immune subtypes, we performed unsupervised consensus clustering on the 26 overlapping TSGs in a discovery cohort of 35 samples (in-house + GSE204663, lacking clinical annotations). Optimal cluster number was determined as k = 2 ( Supplementary Figures S1A–C ), resulting in two robust immune subtypes: Subtype A (n = 20) and Subtype B (n = 15) ( Figure 2G ). Subtype A was enriched for CCL3, IL2RA, and IL1R2, suggesting a T cell–dominant immune response. In contrast, Subtype B showed upregulation of IRF4 and TNFRSF17, genes involved in B cell activation and terminal plasma cell differentiation ( Supplementary Figure S1D ).

To investigate the functional divergence between TLS-based immune subtypes, we conducted differential expression analysis followed by Gene Ontology (GO) enrichment using ClueGO. The resulting functional network revealed 11 significantly enriched GO terms ( Figures 3A, B ), most of which were related to immune processes.

Importantly, Subtype B demonstrated marked enrichment of B cell-associated pathways, including “B cell differentiation,” “positive regulation of immunoglobulin production,” and “mature B cell differentiation.” These pathways were driven by genes consistently upregulated in Subtype B ( Supplementary Table S1 ), reflecting a transcriptional program indicative of active humoral immunity.

Notably, 36.36% of the top enriched GO terms were directly linked to B cell function ( Figure 3C ), supporting the notion that Subtype B is characterized by a B cell–dominant immune phenotype that may underlie its association with more aggressive clinical behavior.

We next assessed the clinical relevance of TLS-associated subtypes using the GSE213346 dataset, which includes clinical annotations for 40 OLP patients. Patients were categorized as non-erosive OLP (NEOLP, n = 27) or erosive OLP (EOLP, n = 13). Following one year of follow-up, they were further stratified into a recurrent erosion (RE) group defined by erosive episodes recurring within 3 months and a persistent non-erosion (PNE) group, who remained erosion-free during the same period, as previously described (3). For comparison, 10 healthy controls from our in-house RNA-seq dataset were included.

Subtype B-associated genes IRF4 and TNFRSF17 were significantly upregulated in both EOLP and RE groups. Similarly, canonical B and plasma cell markers CD79A, CD79B, CD19, MZB1, and CD38 were also elevated in patients with erosive and recurrent disease ( Figures 4A, B ). In contrast, among Subtype A-enriched genes, IL2RA was associated with both erosion and recurrence, whereas CCL3 was linked only to erosion. IL1R2 and CD8⁺ T cell markers (CD8A, CD8B) showed no significant correlation with clinical outcomes ( Figures 4A, B ). These findings reinforce the role of humoral immune activation-particularly B cell maturation and plasma cell differentiation-in driving aggressive OLP phenotypes characterized by erosion and recurrence.

To validate transcriptomic findings at the protein level, we performed immunohistochemical analysis of CD20 (B cell marker) and CD38 (plasma cell marker) in an independent retrospective clinical cohort of 74 patients (26 OLP, 48 OLL).

CD20⁺ B cells were detected in 97.3% of cases and were significantly associated with TLS positivity (P < 0.001), but showed no correlation with clinical subtype or severity ( Table 1 ; Figure 5 ).

In contrast, CD38⁺ plasma cells were present in 47.3% of cases and demonstrated a significant association with erosive disease (P = 0.038). Higher grades of plasma cell infiltration were more frequently observed in erosive rather than reticular lesions, and REU scores tended to increase with CD38⁺ cell density (P = 0.065). Notably, OLL cases were overrepresented in higher plasma cell infiltration grades (P < 0.001), while TLS positivity remained uniformly high across all plasma cell levels (P = 0.510) ( Table 2 ; Figure 6 ).

These results suggest distinct roles for B cell subsets in OLP/OLL pathogenesis: CD20⁺ B cells contribute to TLS formation and structural organization, whereas CD38⁺ plasma cells reflect terminal humoral activation and are associated with more severe clinical phenotypes.