Stereo-XCR-seq achieves high-fidelity spatial immune repertoire mapping with robust performance in fresh frozen and FFPE disease tissues
To acquire an unbiased spatial immune repertoire in parallel with subcellular-resolution transcriptomes, we developed Stereo-XCR-seq, a method that integrates Stereo-seq16 with targeted retrieval of full-length V(D)J sequences (Fig. 1a, Methods). Firstly, double-stranded Stereo-seq cDNA was heat-denatured at 95 °C and rapidly annealed to a splint oligonucleotide, yielding a single-stranded circular intermediate (sscirDNA) in which the barcode resides at the 5′-end (Fig. 1a, step 1–2). We then sealed the nick using T4 DNA ligase, and typically 10 ~ 30% of the input cDNA was converted to circles. Residual linear DNA and excess oligos were removed by exonucleases I and III, eliminating non-specific amplification in subsequent steps (Fig. 1a, step 2). XCR transcripts were then amplified from the sscirDNA pool by PCR primers directed to the constant (C) region (sscirPCR; Fig. 1a, step 3). In contrast to multiplex V-gene PCR, this constant-region strategy is repertoire-agnostic and yields unbiased amplicons in which barcodes and V(D)J segments are in the middle and flanked by truncated C regions. These amplicons were re-circularized for paired-end 150 bp (PE150) sequencing: read 1 captures the CDR3, while read 2 contains the CID and unique molecular identifier (Fig. 1a, step 4). The paired sequencing data were then subjected to barcode mapping and clone assembly, respectively, to examine the spatial distribution of each clone read (Fig. 1a, step 5). Notably, we also provided an optional fragmentation pipeline to substitute step 4 (Supplementary Fig. 1a, methods). The fragmented V(D)J reads sharing identical coordinates and unique molecular identifier could be further assembled to generate full-length V(D)J reads (Supplementary Fig. 1a).
Fig. 1: Stereo-XCR-seq achieves high-fidelity spatial immune repertoire mapping with robust performance in fresh frozen and FFPE disease tissues.
a Schematic graph shows five key steps of Stereo-XCR-seq. Created in BioRender. Zhan, X. (2026) https://BioRender.com/9bbcvwx. b Paired box-strip plot comparing pre-retrieval (Stereo-seq) and post-retrieval (Stereo-XCR-seq) target reads/total reads ratios in fresh frozen tissues. Box-and-whisker plots display median, interquartile range and minimum to maximum values. Each dot represents an independent experimental replicate. The number of replicates (N) is labeled on the graph, n = 8 (TRAC), 9 (TRBC),7 (TRDC), 7 (TRGC), 7 (IGHA/M/D/E/G), 11 (IGKC), 13 (IGLC) samples. Lognormal data were analyzed by two-tailed paired Student’s t test, *p < 0.05, **p < 0.01, ****p < 0.0001. c Paired strip plots show ratio differences among three enrichment approaches. Connected lines represent identical cDNA libraries, N = 3 technical replicates. Two-tailed paired Student’s t test showed significance between sscirPCR and the other two strategies, *p < 0.05, **p < 0.01. d Bar chart shows the CDR3 clone types of each tissue section. Data are presented as mean ± SEM, lognormal data were analyzed by two-tailed unpaired Student’s t tests, *p < 0.05, **p < 0.01, ****p < 0.0001. Biological replicates (N) are labeled in the figure, n = 9 (Stereo-XCR-seq), 5 (Spatial-VDJ), 4 (SPTCR-seq), 4 (Slide-TCR-seq), 1(Slide-tags) samples. e. Bar chart shows the TCRβ and IgH chain cell number of each tissue section using stereo-XCR-seq. Data are presented as mean ± SEM. Biological replicates (N) are labeled in the figure. N = 9 samples. f Representative spatial plots show the density and distribution of XCR transcriptomic reads and clone reads at a resolution of bin50, scale bar = 1 mm. N analyzed samples = 11. For data presentation, IGHA/D/E/M/G are summed as IGH gene and TRBC1/2 are summed as TRBC gene. Clone reads are presented using color map Oranges, whereas transcriptomic reads are presented using color map Greens. Both colors are shaded by UMI counts. FOV selected is highlighted in white frame, scale bar = 100 μm. g Representative mIF staining of CD3E and CD138 using adjacent 5 μm tissue section. N analyzed samples = 3. FOVs corresponding to (f) are shown, scale bar = 100 μm. h Representative spatial profiling of XCR gene expression and clonal reads at bin20 resolution, scale bar = 1 mm. TCR and BCR signals are distinguished by color and shaded by CID and UMI counts. N analyzed samples = 3. i Paired box-strip plots of read ratios in FFPE samples. Box-and-whisker plots display median, interquartile range and minimum to maximum values. Each dot represents an independent experimental replicate. N = 5 samples. Lognormal data were analyzed by two-tailed paired Student’s t test, **p < 0.01, ***p < 0.001, ****p < 0.0001. Source data are provided as a Source Data file.
Although long-read platforms recover complete V(D)J-C sequences, their base-level accuracy is limited. Conversely, short reads provide high fidelity but only partial V(D)J coverage. Previous studies have utilized both sequencing methods to mitigate the shortcomings of each method18, whereas the readouts of the two methods have not been compared in parallel. Thus, we directly compared nanopore-based long-read sequencing with DNA nanoball-based PE150 short reads for spatial XCR-seq. Long reads exhibited significantly lower quality (Supplementary Fig. 1b) and produced 39.9% non-functional T/BCR sequences, versus 6.9% for short reads (Supplementary Fig. 1c). In contrast to the plateaued curve of short-read sequencing, increasing long-read sequencing throughput led to ever-increasing CDR3 clone types per coordinate in the saturation test, indicating accumulated false positive readouts by sequencing errors (Supplementary Fig. 1d). These results strongly suggested that short-read sequencing result is imperative to generate high-fidelity clonotype information in spatial immune repertoire. To obtain full-length information, we reconstructed full-length V(D)J through the fragmentation pipeline (Supplementary Fig. 1a). Coverage analysis showed a high coverage ( > 90%) of CDR1-to-frame region 4 (FR4) in the assembled clone reads (Supplementary Fig. 1e). Compared to previous studies that utilized the advantages of both sequencing methods18, our method provides a short-read-only strategy to generate full-length clone reads that enhance the fidelity of the sequencing results.
To examine the performance of Stereo-XCR-seq, we applied sscirPCR to multiple Stereo-seq cDNA libraries and retrieved XCR reads with consistent efficiencies (5.2 ~ 13.5 log2 fold changes over the original Stereo-seq cDNA library, Fig. 1b). To demonstrate the robustness of sscirPCR, we benchmarked the clones detection against multiplexed PCR and 5’RACE19 (rapid amplification of cDNA 5’ ends, a widely recognized gold standard for bulk and single‑cell TCR/BCR sequencing, Supplementary Fig. 1f). Compared to previously released spatial T/BCR retrieving strategies, such as multiplexed PCR12,20 and probe hybridization15, sscirPCR showed substantially higher retrieving efficiencies in XCR reads enrichment (Fig. 1c). The significantly correlated clone detection by sscirPCR and 5’RACE on the same tissue section confirmed the consistency with the gold‑standard approach (r = 0.881, p = 6 × 10−8, Supplementary Fig. 1g). Given the extreme diversity of human T/BCR repertoires, reliance on predefined V(D)J primer panels can introduce substantial coverage bias across samples. In such a case, targeting the constant region might be a more optimal solution to reduce this bias. To prove the unbiased enrichment by sscirPCR, we conducted sscirPCR, 5’RACE, and multiplexed PCR to enrich the TCR reads using the same cDNA library for orthogonal validation. The pairwise clone detection correlation analysis showed significant linear correlation between 5’RACE and sscirPCR (Supplementary Fig. 1h, left), while the multiplexed PCR exhibited marginal correlation with either 5’RACE or sscirPCR (Supplementary Fig. 1h, middle and right). In addition, sscirPCR and 5’RACE have very similar TRBV gene coverage (Supplementary Fig. 1i, left). Multiplexed PCR detected 22 V genes, nearly all of which were also identified by the other tw. methods (Supplementary Fig. 1i, middle and right). This comparison suggests that sscirPCR could achieve unbiased enrichment.
Owing to its high efficiency, sscirPCR was also able to retrieve TCRγ and δ chains in addition to TCRα and β chains, even though γδT cell only accounts for a small fraction ( ~ 5%) of CD3+ T cells21,22. In addition, we tested sscirPCR across a broad range of tissue types, including LUAD, colorectal, gastric, bladder, renal, lymph nodes with metastatic esophageal cancer, autoimmune hepatitis, and rheumatoid arthritis. On average, we assembled 2145 TCRβ and 4659 IgH clones that cover 10,845 TCRβ+ and 51,639 IgH+cells per sample library (Fig. 1d, e), demonstrating robust performance of sscirPCR. Compared with published spatial immune repertoire datasets, Stereo-XCR-seq exhibited a strong advantage in the CDR3 clone type counts23 (Fig. 1d). Furthermore, the retrieved XCR clone reads exhibited consistent spatial distribution along with original spatial transcriptomic T/BCR gene expression (Fig. 1f). The existence of T cells and plasma cells in highlighted T/BCR-expressing field of views (FOVs) was further validated by CD3E and CD138 immunofluorescence staining (Fig. 1g). Importantly, Stereo-XCR-seq also displayed capabilities in XCR sequencing using FFPE samples from colorectal cancer and breast cancer (Fig. 1h, i, source data file), though notable background noises were observed before filtering (Supplementary Fig. 1j). Due to the harsh preservation condition, the XCR reads in FFPE tumor specimens from breast cancer were degraded severely, leading to a shorter median XCR read length at 279 bp (versus 1329 bp for fresh frozen tissues) that only covered the CDR3-FR4 region (Supplementary Fig. 1k–l). Nevertheless, Stereo-XCR-seq still obtained 1026 ~ 2840 IgH and 240 ~ 593 TCRβ clones per sample (source data file), thus extending the application in clinical studies using retrospective FFPE samples. Collectively, the high-fidelity readout, enhanced retrieval efficiencies, unbiased enrichments, and robust performances in different tissues under varied preservation conditions highlight the application potential of Stereo-XCR-seq in exploring lymphocyte dynamics in immune-related studies.
Stereo-XCR-seq profiles clone activities of lymphocytes at the single-cell level with paired chains
Effective anti-tumor immunity requires clonal expansion of tumor-infiltrating lymphocytes (TILs)24, accurate quantification of clone sizes is therefore essential for inferring tumor reactivity. To achieve single-cell resolution, we incorporated ssDNA staining into the Stereo-seq workflow, capturing high-resolution nucleus staining images along with spatial transcriptomes (Supplementary Fig. 2a). Using cellbin2 software25, we performed image-based cell segmentation (Supplementary Fig. 2b). By manually checking 30 FOVs randomly selected from three LUAD samples (10 FOVs each), we obtained authentic segmentation rates at ~94%, ambiguous segmentation rates at ~6%, and negligible missing cells (Supplementary Fig. 2c, d), supporting a credible cell segmentation through this pipeline. A critical consideration in cell segmentation is avoiding excessive boundary expansion, which risks incorporating spurious signals from neighboring cells. To examine this, we first quantified the endogenous membrane-to-nucleus area ratio of T lymphocytes. On average, the membrane area was 1.687-fold larger than the nuclear area (n = 30 cells, Supplementary Fig. 2e, f, source data). For T-cell segmentation, our expanded cell masks were 1.684-fold larger than the original nuclear masks—a fold change nearly identical to the membrane-to-nucleus ratio, with no statistically significant difference (Supplementary Fig. 2f). This result indicates that the 10 pixel expansion parameter does not lead to excessive boundary expansion. Existing spatial immune-repertoire methods typically have a resolution at 100 µm (center-to-center distances, such as Spatial VDJ, SPTCR-seq)26, while spots at this resolution could contain a varied number of lymphocytes (Fig. 2a). Quantifying TCRβ⁺ and IgH⁺ cells revealed a mean of 2.58 TCRβ⁺ and 6.58 IgH⁺ cells per 100 µm spot (Fig. 2b), suggesting an inaccurate clone size counting at a coarse-grained resolution. We randomly selected 20 pairs of FOVs from three samples for the Stereo-XCR-seq section and the serial section stained by CD3E/CD20/CD138. The quantification results of 100 µm XCR+ spots and XCR+ cell bins were directly compared with the mIF staining ground truth in the serial section. As a result, the cell bin counts of either T cells or B/plasma cells well correlated with the cell counts of the mIF staining results, whereas cell counts by 100 µm (bin 200) spots showed no significant correlation (Fig. 2c, left). We then plotted T-cell and plasma cell counts of mIF staining and of Stereo-XCR-seq at 100 µm resolution and cell bin resolution in three pairs of FOVs (Fig. 2c, right). In contrast to the deviated quantification of 100 µm spot counts, cell bin exhibited well-correlated T-cell and plasma cell counts with the mIF staining cell counts on the serial tissue sections. Collectively, the above result suggests that Stereo-XCR-seq provides credible clone size quantification at single cell resolution.
Fig. 2: Stereo-XCR-seq profiles clone activities lymphocytes at single-cell level with paired chains.
a Spatial plot shows the distribution of IgH clone family 368 at two different resolutions, 100 μm and single-cell resolution. FOV is highlighted in white frame. Bin200 spot (100 μm) is plotted as pink frame. Single cells are plotted in pink polygons with white cell border. b Violin-box plots show the TCRβ+ and IgH+ cell counts in each bin200 spot. Box-and-whisker plots display median, interquartile range and minimum to maximum values. Biological replicates (n) represent different spots from the same tissue section. c Two-tailed Pearson’s correlation analysis (left and middle) using T-cell and B/plasma counts at cell bin resolution, bin200 resolution (100 μm) and mIF staining. N = 3 samples. Each dot represents an individual FOV (N = 20 FOVs), with colors indicate different samples. r- and p values are labeled on each analysis. T-cell and B/plasma cell counts from three FOVs are plotted to show the counts distribution of cell bin resolution, bin200 resolution (100 μm) and mIF staining (right). Dots are colored by sample and shaped by data type. d Schematic diagram shows the definition of ambiguous paired chains, unpaired chains, and unique paired chains. e Pie charts show the number of cells/bin200 spots with ambiguous paired chains, unpaired chains, and unique paired chains. f, g Pie chart (f) and Spatial plot (g) shows the capacity of Stereo-XCR-seq in resolving clone expansion. h Spatial plot shows the capacity of Stereo-XCR-seq in resolving CSR events. each dot represents a transcript, colored by isoforms. i Getis-Ord Gi* hotspot map (k = 10) showing statistically significant spatial aggregation of cells bearing highly similar IgH CDR3 sequences. Insets highlight examples of significant spatial clustering. j Lineage tree plot and spatial plot shows the clone family 340 and subclones 2/6/8 in different FOVs. Only subclones 2/6/8 are shown in spatial plot. In lineage tree, dots are sized by clone size. In spatial plot, each dot/polygon represent a cell, colored by subclone definition. Source data are provided as a Source Data file.
Equipping T cells with chimeric antigen receptor (CAR) or screened TCR endows them with tumor-targeting cytotoxicity in clinical practice27,28. However, such engineering requires information on paired TCR chains. To this end, we compared the resolution of paired chain detection at single cell level and at 100 µm. By examining the co-localization of TCR α-β and BCR H-L chains within cell bins and 100 µm spots, we categorized the result into ambiguous paired chains with multiple distinct clonotypes, unpaired chains with only one chain type, and unique paired chains with a single, clearly defined clonotype (Fig. 2d). Compared to the 100 µm resolution, Stereo-XCR-seq exhibited no ambiguous pairings and yielded a higher proportion of unique paired chains at single cell resolution (n = 8785, 61.6%, Fig. 2e). Collectively, these result demonstrate that Stereo-XCR-seq enhables single-cell spatial immune repertoire profiling with paired chains, offering potentials for retrieving paired T/BCR receptors from clinical specimens to support cellular therapies.
B/plasma lymphocytes undergo somatic hypermutation (SHM) to develop antigen-binding affinity29, class switching recombination (CSR) to secrete effector-competent antibodies, and clonal expansion to enable robust antibody secretion30. These three clonal activities are essential for B/plasma cells to fulfill their roles in immune surveillance. We therefore sought to evaluate the capability of Stereo-XCR-seq in resolving each of these clonal activities at single-cell resolution. To assess clonal expansion, we quantified the number of cells assigned to each clone family to determine clone sizes (Fig. 2f). For example, in LUAD-P1, 11% (n = 71) IgH clones were found with large clone size ( ≥ 10 cells, such as clone family 425 in Fig. 2g). For CSR, we defined a CSR event as the presence of two or more IgH isotypes within a segmented cell bin (Fig. 2h). Using sequence similarity analysis and Getis-Ord Gi* statistics, we identified a mutation hotspot characterized by similar IgH CDR3 amino acid sequences in this FOV (Fig. 2i, Supplementary Fig. 2g). By integrating these capacities, Stereo-XCR-seq enables the visualization of individual clonal activities involved in B/plasma-cell lineage differentiation (Fig. 2j, left, orange bold lines). For example, in the IgH clone family 340, clonal expansion of subclone 2 (clone size at 77), mutation from subclone 2 to subclone 8, and a CSR event from subclone 2 to subclone 6 could all be inferred within a single spatial plot (Fig. 2j, right), demonstrating the capability of Stereo-XCR-seq to resolve B/plasma cells clonal dynamics at single cell resolution.
Stereo-XCR-seq reveals clonal convergence of B/plasma cells across spatially distinct iTLSs in LUAD
TLSs are ectopic lymphoid aggregates within tumor stroma that orchestrate anti-tumor immunity31,32. The maturation of TLS is marked by the formation of germinal centers (GCs) that generate antibody-secreting plasma cells within the follicular structures33,34. The presence of mature TLS (mTLS) has been reported to associate with improved response to immunotherapy. However, only a minority of LUAD patients harbor mTLSs35,36,37. To investigate their prevalence, we randomly selected H&E staining images from 100 LUAD samples in The Cancer Imaging Archive38 (Supplementary data 1). Of these, 57 samples contained five or more TLSs, indicative of an inflamed tumor microenvironment. However, only six samples displayed mTLSs (Fig. 3a, b). The mean frequency of mTLS remained significantly lower than that of iTLS across the inflamed cohort (mean frequencies: 0.26 vs 10.05, Fig. 3c), substantially below the reported objective response rates to immune checkpoint blockade in LUAD (35.9%39 and 30%40). These findings suggest that the prognostic value of mTLS applies only to a subset of LUAD patients.
Fig. 3: Stereo-XCR-seq reveals clonal convergence of B/plasma cells across spatially distinct iTLSs in LUAD.
a Representative FOVs show the iTLSs and mTLSs. Germinal centers are highlighted in green. Patient number has been attached below each image. N analyzed samples = 100, scale bar = 100 μm. b Pie chart shows the percentage of mTLS and iTLS appearance in LUAD H&E staining images. c Box plot shows the frequencies of mTLS and iTLS in 57 inflamed LUAD H&E staining images. Box-and-whisker plots display median, interquartile range and minimum to maximum values. Statistical analysis was performed using two-tailed unpaired Student’s t tests. ****p < 0.0001. d Schematic diagram shows the investigation of LUAD-infiltrating lymphocytes at single-cell resolution and structural resolution. Created in BioRender. Zhan, X. (2026) https:// BioRender.com/9bbcvwx. e Lollipop charts show the number of TCR+ cells and BCR+ cells in the 11 LUAD tumors. f Representative spatial plot shows iTLS in LUAD-P1, n = 7338 cells, scale bar = 1 mm. FOVs presents the region of interests, along with the CXCL13, MS4A1, CD3E expression in this FOV. Each dot represents a bin50 spot, colored by cluster or expression level respectively. N analyzed samples = 11. g Spatial plot shows spatially discrete iTLS clusters in LUAD-P1, n = 6243 cells, scale bar = 1 mm. N = 6 samples (LUAD-P1/P4/P6/P8/P9/P10, with iTLS number ≥2). h Venn plots show the shared IgH clones by the top five iTLS clusters in LUAD-P1. iTLSs 7/3/8/0 are compared with iTLS 6, respectively. i Elbow plots show the sharing of IgH clones (left) and TCRβ clones (right) by different iTLS clusters, n = 6 samples (LUAD-P1/P4/P6/P8/P9/P10, with iTLS number ≥2). The sharing clones are plotted as absolute counts (up) and percentage (below), respectively. Each line represents a tumor, colored by patient number as indicated on the left. j Box-whisker chart (left) and schematic graph (right) shows the disseminated clone proportions in each iTLS. Box-and-whisker plots display median, interquartile range and minimum to maximum values. Biological replicates (N) are labeled in the figure, n = 10 samples (LUAD-P1 to P10, with iTLS number ≥ 1). k Representative spatial plots show the distribution of iTLS-emigrating IgH clones, scale bar = 1 mm. Each plot presents one iTLS cluster-associated clones as labeled above. The iTLS cluster of interest is colored in blue and highlighted by red border. N analyzed samples = 6 (LUAD-P1/P4/P6/P8/P9/P10, with iTLS number ≥2). The emigrating distributions are presented as kernel density estimate (KDE) plot, with the most enriched niche highlighted by green borders. The KDE plots are colored by densities. FOVs show the representative distribution of cells of interests and related gene expression in LUAD-P1. Each dot represents a cell, colored by either cell type or genes. Dot plot shows the expression of JCHAIN in different spatial clusters (bin50 spots). Dots are colored by mean expression and sized by expression fraction. Source data are provided as a Source Data file.
To investigate the role of iTLSs in LUAD, we applied Stereo-XCR-seq to 11 fresh-frozen LUAD specimens containing only iTLSs, aiming to characterize the clonal activity of lymphocytes within spatially distinct lymphoid aggregates (Fig. 3d, Supplementary Table 1). Serial tissue sections were used for H&E staining and LAMP3/CD23 co-immunostaining. Although LAMP3 expression was observed in the TLSs, the absence of CD23+ secondary follicular structures confirmed their immature phenotype41 (Supplementary Fig. 3a–c). Unsupervised clustering analysis42 (built in Omicverse43) of spatial transcriptome at 25 μm × 25 μm (bin50 spot) resolution resolved seven spatial structures (Fig. 3f, Supplementary Fig. 3d, e): alveolar tissue (SFTPC+ SFTPD+), lung cilia (SCGB1A1+CAPS+), plasma cell aggregates (JCHAIN+IGKC+), iTLS (CXCL13+MS4A1+CD3E+CHST4+IL33+), myeloid cell-enriched region (APOC1+ CTSB+), immune-excluded stroma (CXCL13−MS4A1−JCHAIN−COL1A1+COL1A2+) and tumor (EPCAM+KRT18+). In addition, we transferred cell-type annotations from a public LUAD scRNA-seq atlas (GSE148071)44 to our in-house generated single-cell resolved spatial data using TACCO45. Marker-gene concordance between scRNA-seq data and spatial transcriptomic cell bins validated the accuracy of the annotation transfer (Supplementary Fig. 3f–h). Leveraging the single-cell resolved transcriptome and spatial annotations at 25 μm resolution, we were able to link individual gene expression profiles to their respective tissue niches.
Given that single-cell resolved spatial transcriptome is challenged by background noises from lateral diffusion and the Z-axis interference46, we restricted our analysis to XCR reads assigned to lymphocytes to ensure biological validity in downstream analysis (Supplementary Fig. 4a, b). Specifically, we recovered 14,125 TCR+ T-cell (recovering rate at 66.9%) and 22,084 BCR+ B/plasma cells (recovering rate at 73.33%) (Fig. 3e), encompassing 520 TCRβ clones and 3130 IgH clones. Among these, ~88% (17,946/20,361) IgH+ cells were supported by IGH gene expression, whereas only 9.7% (757/7,826) of TCRβ+ cells exhibited detectable TRBC1/2 gene expression (Supplementary Fig. 4c–e). A recent study comparing Xenium and Visium platforms demonstrated that sequencing-based spatial transcriptome based whole-transcriptome methods underestimates the expression of low-abundance immune-related genes by 5.7-fold per gene17. We thus attributed the discrepancy in supporting ratios (Supplementary Fig. 4e) to the relatively low abundance of TCR mRNAs in spatial transcriptome cDNA libraries. This interpretation is further supported by our own data and other publicly available LUAD spatial trnacriptomes (E-MTAB-13526)47, which consistently showed lower TRBC1/2 mRNA levels compared to IGH gene expressions (Supplementary Fig. 4f, g). This observation prompted us to question whether TCRβ+ cells lacking TRBC1/2 expression represented biologically authentic cells restored via retrieval or false positive signals introduced from the sscirPCR procedure. To test this, we checked whether TCRβ+TRBC1/2− cells were supported by TRAC mRNA, TCRα chain reads, or CD3E protein expression. Although lacking TRBC1/2 expression, TCRβ+TRBC1/2− cells showed higher TRAC expression compared to TCRβ-TRBC1/2− cells (Supplementary Fig. 4h). Notably, ~67% TCRβ+ cells harbored paired TCRα chains, resulting in 5202 T cells with recovered α-β chains (Supplementary Fig. 4i). We further assessed CD3E expression at protein level by immunofluorescent staining using serial tissue sections. We selected two representative fields of view (FOV), including FOV1 with concordant TRBC1/2 and TCRβ expression, and FOV2 with discordant expression between the two (Supplementary Fig. 4j). Comparable CD3E expression in both FOVs (Supplementary Fig. 4k) supported the authenticity of the retrieved TCRβ signals, highlighting the enhanced sensitivity of Stereo-XCR-seq. Overall, the T/BCR+ cell counts showed strong agreement with the number of T/B lymphocytes per tissue section (Supplementary Fig. 4l) and at 100 μm2 regional resolution (Supplementary Fig. 4m), underscoring the accuracy of clonal recovery.
Among these specimens, six tumors (LUAD_P1/4/6/8/9/10) contained multiple iTLSs. To identify spatially distinct iTLS clusters in these tumors, we applied Density-Based Spatial Clustering of Applications with Noise (DBSCAN), a density-based algorithm designed to detect clusters in large spatial databases with noise48. As a result, we identified 3 ~ 10 spatially discrete iTLS clusters per tumor section (Fig. 3g, Supplementary Fig. 5a). By checking the shared IgH clones by the top five iTLSs in LUAD_P1, we observed a minor clone sharing ranging from 11.7 ~ 27.3% (shared clones/all clones, Fig. 3h), implicating spatial heterogeneity of iTLSs. To further validate this observation, we calculated the clone type sharing by different iTLSs in all six LUAD tumors with multiple iTLSs and plotted the result as elbow plots (Fig. 3i). Again, ~70% IgH clones and ~60% TCRβ clones were found to be unique (appeared in only one iTLS) in different iTLSs of each tumor, whereas <5% IgH clones and <10% TCRβ clones were found shared by three or more iTLSs within each tumor (Fig. 3i). Of note, the pairwise clone sharing analysis of the iTLSs showed no significant correlation between clone sharing decay with distance, indicating that clonal overlap might not be directly determined by spatial proximity (Supplementary Fig. 5b). Interestingly, though rare clones were shared by spatially discrete iTLSs, we observed that over 80.67% IgH and 60.61% TCRβ clones could also be observed outside the iTLSs, suggesting a crucial role of iTLS in supplying tumor-infiltrating lymphocytes (Fig. 3j). Next, we integrated the transcriptome and clone information to interpret the migration paths of B/plasma cells in the tumor. Mapping iTLS-associated IgH clones revealed extrafollicular egress and convergence on discrete stromal niches (Fig. 3k, dotted arrows), which was further enhanced by the distant aggregation of iTLS-associated IgH clones in all examined multi-iTLS tumors (Supplementary Fig. 5c). B cells differentiate into plasma cells to secrete anti-tumor antibody23. While iTLSs harbored dense B and T cells, the emergent niches were overwhelmingly occupied by plasma cells, as evidenced by the intensive expression of JCHAIN in the distinct niches (Fig. 3k, left FOVs and dot plot). This observation indicates that the emigrant B cells complete plasma cell differentiation at these sites. Given the aggregation of plasma cells in the distinct niche, we named it plasma cell zones (PCZ). Collectively, our data suggested that the fate transition from B cells to antibody-secreting plasma cells might be a process that proceeds asynchronously during migration.
Stereo-XCR-seq delineates spatiotemporal fate trajectories of B/plasma cells in iTLSs
To resolve whether B-cell maturation is a spatially dynamic process, we projected MS4A1, CXCL13, IGKC, and JCHAIN onto the spatial maps from the 11 LUAD tumors. It turned out that iTLSs and PCZs showed distinct spatial distribution patterns marked by co-expression of MS4A1/CXCL13 or IGKC/JCHAIN, respectively (Fig. 4a, b). Cell-type deconvolution and CD20/CD138 co-staining on serial tumor sections further confirmed the enrichment of B cells in iTLSs and plasma cells in PCZs, respectively (Fig. 4c, d, Supplementary Fig. 6d–h). Combining the observation of shared IgH clones in iTLSs and PCZs (Fig. 3k), the distinct intratumoral distributions of B cells and plasma cells implied that B-to-plasma cell differentiation might be coupled to intratumoral relocation. By generating PCZ and TLS gene signature scores (Methods), we found that this observation extended to immune disorders and cancers in other tissues, including primary biliary cholangitis (PBC), inflammatory bowel disease (IBD), breast cancer (BC), colorectal cancer (CRC), kidney cancer (KC) and gastric cancer (GC). These conditions exhibited distinct spatial expression patterns of TLS and PCZ scores, along with related gene markers (Supplementary Fig. 6a–c). The conserved spatial distribution of iTLSs and PCZs across diverse tissues and diseases not only supports the biological relevance of our findings in LUAD but also indicates that B-to-plasma cell maturation could occur outside iTLS at a pan-disease scale.
Fig. 4: Stereo-XCR-seq delineates spatiotemporal fate trajectories of B/plasma cells in iTLSs.
a Spatial plot of MS4A1 (green), CXCL13 (red), IGKC (purple) and JCHAIN (blue) expression in LUAD tumor, bin20 spot, colored by genes and UMI counts, n = 11 samples, scale bar = 1 mm. Two FOVs are selected to present the co-expression of MS4A1/CXCL13 (green frame) or IGKC/JCHAIN (purple frame). b Spatial plot of selected FOVs from (a) showing co-expression of MS4A1/CXCL13 (green frame) or IGKC/JCHAIN (purple frame) in bin20 spot, colored by genes and UMI counts. n = 11 samples, scale bar = 10 μm. c, d Representative multiplexed immunofluorescence staining (c) of LUAD-P6 on an adjacent tissue section shows the CD20+ B cells and CD138+ plasma cells, scale bar = 1 mm. FOVs are selected according to the location of PCZ and iTLS to show the B cell enrichment in iTLS and plasma cell enrichment in PCZ, scale bar = 20 μm. Box-strip plots (d) compare their abundance; two-tailed paired Student’s t tests, **p value < 0.01. Box-and-whisker plots display median, interquartile range and minimum to maximum values. N = 3 samples, n = 5 FOVs of each sample. e Schematic of niche/cell analysis at bin50 and single-cell resolution. f UMAP of lymphoid aggregate re-clustering (left-up) and gene expression (other plots, bin50 spots, colored by cluster/UMI), n = 11 samples, n = 16398 cells. g, h Representative spatial plot of inner and outer iTLSs in LUAD P4 (n = 3941 cells) and P6 (n = 4815 cells), bin50 spot, colored by clusters (g) or gene UMI counts (h). Scale bar = 1 mm. N analyzed samples = 11. i Dot plot of gene expression in iTLSs and PCZs (bin50 spots), n = 11 samples. j, k UMAP with slingshot pseudotime trajectory (j) and minimum spanning tree plot (k) shows the B-to-plasma cell fate transition, n = 11 samples, n = 21,682 cells. Each dot represents a cell, colored by pseudotime score (j) and spatial distribution (k). l Representative spatial quiver maps show the migrating paths of B/plasma cells in iTLS clusters (LUAD-P4, n = 233 cells; LUAD-P6, n = 535 cells). B/plasma cells are shown as dots colored by spatial cluster. Black vectors represent the migrating path; scale bar = 1 mm. FOVs are selected to present the B/plasma cells egress. N analyzed samples = 11. m Representative spatial single-cell plot shows the pseudotime scores of B/plasma cells in different niche of LUAD-P6 (n = 219 cells), scale bar = 1 mm. Each dot represents a cell, colored by niches. FOV shows the clone family 101 of the framed region at single-cell resolution. Each polygon represents a cell, colored by pseudotime score. N analyzed samples = 11. n Dot plot of LUAD-P6 top 10 clone families’ pseudotime scores in different niches. Source data are provided as a Source Data file.
We next reconstructed the molecular itinerary underpinning this relocation, aiming to elaborate on the spatiotemporal fate transition from B cells to plasma cells in LUAD. To this end, we analyzed B/plasma cells at cell bin resolution and their corresponding niches at 25 μm resolution (bin50 spots), linking cell fate transition to their spatial distribution (Fig. 4e). Firstly, we re-clustered bin 50 spots corresponding to iTLSs and PCZs across 11 LUAD specimens, as ~70% all B/plasma cells were concentrated within these two structures ( ~ 60% in PCZ, ~10% in iTLS, Supplementary Fig. 7a). Unsupervised secondary clustering subdivided these two structures into inner iTLS, outer iTLS, IgM+ PCZ, and IgG+ PCZ (Fig. 4f). IgM+ PCZ and IgG+ PCZ exhibited exclusive expression of IGHM and IGHG1/2/3/4, respectively (Fig. 4f, Supplementary Fig 7b). Inner iTLS and outer iTLS shared expression of MS4A1, while inner iTLS was surrounded by outer iTLS (Fig. 4f, g). As expected, inner iTLSs displayed higher human lymphocyte antigen (HLA) class I (HLA-A/B) and class II (HLA-DQB1/DPA1/DPB1) expression, indicative of active antigen presentation, whereas outer iTLSs exhibited elevated JCHAIN, signifying ongoing plasma-cell commitment (Fig. 4h and Supplementary Fig. 7c). This observation suggested distinct roles of inner iTLSs and outer iTLSs by harboring different activities in immune activation.
Next, we examined B/plasma cells residing in different niches (Fig. 4e). To assess their differentiation states, we plotted the expression of key marker genes PAX5, XBP1, IRF4, PRDM1, and JCHAIN, which represent sequential stages in B-to-plasma cell differentiation process49,50,51,52,53,54. We compared expression patterns of these markers in B/plasma cells located in iTLSs and PCZs. Interestingly, B/plasma cells in the inner iTLSs exhibited the highest expression of PAX5, suggesting an early stage of differentiation characteristic of centrocytes49,50 (Fig. 4i). B/plasma cells displayed a progressive transcriptional gradient of XBP1, IRF4, PRDM1, and JCHAIN extending from the inner iTLSs to the IgG+ PCZ, suggesting that B cells gradually acquired a terminally differentiated plasma cell phenotype during their migration from iTLS to PCZ51,52,53,54 (Fig. 4i). Using slingshot55, we further constructed a B/plasma cell-specific pseudotime using B/plasma cells from iTLS and PCZ and inferred a differentiating trajectory. As a result, cells residing in IgG+ PCZ and inner iTLS exhibited the highest or lowest pseudotime score, respectively (Supplementary Fig. 7d and Fig. 4j, bottom left). By setting the starting point at inner iTLS, we observed two divergent differentiating paths of B/plasma cells from inner iTLS towards IgM+ PCZ or IgG+ PCZ (Fig. 4j). Using Minimum Spanning Tree analysis56, we observed that MS4A1+ B/plasma cells in outer iTLS could further differentiate into IgM+ B/plasma cells or IgG+ B/plasma cells (Fig. 4k, Supplementary Fig. 7b), in line with the observation that both differentiating paths went through the outer iTLSs. To infer migratory directionality of B/plasma cells in iTLS, we calculated quiver fields57,58 based on the pseudotime of each B/plasma cell. Vector flow analysis revealed directional migration from the inner to the outer iTLSs (Fig. 4l), coinciding with elevated B/plasma cell pseudotime at the iTLS margins (Supplementary Fig. 7e). The outward migration of B/plasma cells likely impedes the formation of ectopic GC-like niches within the iTLSs. To further explore whether B-to-plasma cell differentiation is associated with intratumoral relocation, we plotted the top 10 expanded IgH clone families of LUAD P6. As expected, all top 10 IgH clone families exhibited higher pseudotime scores in IgG+ PCZ compared to other locations (Fig. 4m, n), suggesting a spatiotemporal link between B/plasma cell maturation and migration. Taken together, these data revealed that B-to-plasma cell differentiation is coupled with their intratumoral relocation. The emigration of B/plasma cells from iTLS may underlie the impaired formation of GC-like niches within these structures, implicating that terminal maturation of B/plasma cells occurs ectopically outside the iTLS.
Stereo-XCR-seq reveals PCZ as an ectopic GC-like niche for tumor-reactive lymphocyte priming in LUAD
SHM and CSR events are processes typically occurring in GCs, where antigen-specific B/plasma clones are selected to produce antibodies of distinct isotypes with diverse effector functions54,59,60. Though GCs were not found in iTLSs in all our tumor sections (Supplementary Fig. 3a–c), we observed distal disseminations of the later-staged plasma cells and the early-staged B cells from the same clone families (Figs. 3k and 4i). This observation indicated that B cells acquired the antibody-secreting capacity when maturing in the distinct niches. Using TradeSeq61, we obtained 47 pseudotime-associated genes (Fig. 5a). The top-ranked genes and their non-linear fitting curves indicated a strong association between terminal differentiation and IgG secretion (Fig. 5a, b). To uncover where and how B/plasma cells established antibody-secreting capacity stepwise in the GC-missing tumors, we first defined hypermutated B/plasma cells by IgH chain mutation percentage over 5% based on a per-nucleotide basis, referring to a recently published B cell study62 and class-switched B/plasma cells by containing multiple isoforms (Fig. 2h). Spatial plots of hypermutated B/plasma cells, CSR events, and spatial clusters exhibited co-localizations of SHM and CSR activity with PCZs, but not with iTLSs (Fig. 5c–e). The proportion of hypermutated B/plasma cells (Fig. 5f) and class-switched B/plasma cells (Fig. 5g) was both significantly higher in PCZs compared to iTLSs. In addition, we calculated the CSR events between each pair of isoforms and normalized them with cell counts in each tumor. Compared to iTLSs, PCZ exhibited increased normalized CSR frequencies predominantly featuring IgM-IgG and IgA-IgG transitions (Fig. 5h), as presented by spatial plots (Supplementary Fig. 8a) and box-strip plots (Supplementary Fig. 8b). Furthermore, the highest IGHG1/3/4 expression in PCZs over all other spatial clusters (Fig. 5i) suggested that CSR occurs primarily within PCZ to generate IgG-secreting plasma cells. Collectively, the co-occurrence of SHM and CSR implied that PCZ might play a role similar to canonical GCs, which should have occurred in the TLSs.
Fig. 5: Stereo-XCR-seq reveals PCZ as an ectopic GC-like niche for tumor-reactive lymphocyte priming in LUAD.
a Dot plot of Wald’s test result between gene expression and pseudotime scores. b Multiple curve plots show the nonlinear quartic fitting of normalized unique molecular identifier (UMI) counts of each highlighted genes and pseudotime score in B/plasma cells. c, d Representative spatial KDE plots show the SHM and CSR activity densities in LUAD tumors (LUAD-P10, n = 13140 cells; LUAD-P7, n = 8724 cells, n analyzed samples = 11), scale bar = 1 mm. e Representative spatial plot shows the cells in PCZ and iTLS. FOVs are selected to show the distribution of PCZ and iTLS. Each dot represents a cell, colored by spatial clusters. N analyzed samples = 11, scale bar = 1 mm. f, g Paired Box-strip plot of SHM and CSR B/plasma cell proportion in iTLS and PCZ, n = 11 samples. Statistical analysis was performed using two-tailed paired Student’s t tests. Box-and-whisker plots display median, interquartile range and minimum to maximum values. *p < 0.05. h Heatmap shows the normalized CSR frequencies in iTLS and PCZ, n = 11 samples. i Dot plot shows the expression of the gene of interest in different spatial clusters (bin50 spots), n = 11 samples. j Dot and matrix plots show the mutation percentage and mean expression of T–B-cell interaction candidate molecules in different niches, n = 11 samples. k Paired box-strip plots show the comparison of T-cell infiltration in IgG+ and IgM+ PCZs, n = 11 samples. Statistical analysis was performed using two-tailed paired Student’s t tests. Box-and-whisker plots display median, interquartile range and minimum to maximum values. **p < 0.01. l Representative spatial plot (above) shows the cells in PCZ and iTLS (LUAD-P6, n = 2029 cells), n = 11 samples. FOVs are selected to show the distribution of IgG+ PCZ and IgM+ PCZ. Each dot represents a cell, colored by spatial clusters. Representative multiplexed immunofluorescence staining (below) using the adjacent tissue sections shows the CD20+, CD3E+ and CD138+ cells, respectively, n = 7 samples. FOVs are selected according to the spatial plots above to show the T cell infiltration in IgG+ PCZ, scale bar = 10 μm. m Paired box-strip plots show the comparison of T-cell infiltration in IgG+ and IgM+ PCZ, n = 7 samples. T cells were counted by mIF staining data using Image J (Version: 1.54p). Statistical analysis was performed using two-tailed paired Student’s t tests. Box-and-whisker plots display median, interquartile range and minimum to maximum values. ***p < 0.001. n–q Two-tailed Pearson’s correlation analysis of T-cell signature versus mean B/plasma cell pseudotime scores (n), CSR event counts (o), IgG secretion (p) and IgH mutation percentage (q). Each dot represents a PCZ, colored and sized by pseudotime scores (n, o) or CSR event counts (p, q), n = 11 samples. Pearson’s correlation coefficient (r) and p values are labeled in the figure. r Representative partial plot shows the whole slide distribution of T-cell clones present (upper panel, n = 1757 cells) or absent (lower panel, n = 45 cells) in IgG+ PCZ of the LUAD-P6 sample, scale bar = 1 mm. N analyzed samples = 11. s Grouped box-whisker plots show the clone sizes of present and absent T-cell clones in eight tumors with both present and absent clones. Statistical analysis using two-tailed Mann–Whitney U test (independent samples), analysis restricted to samples with clone size >1 in present or absent groups. Box-and-whisker plots display median, interquartile range and minimum to maximum values. Exact p value are provided in source data file. **p value < 0.01, ****p value < 0.0001. t Dot plot shows the expression of gene of interest in the ectopic light zones. u Graphic abstract shows the ectopic maturation of B/plasma cells in LUAD. Source data are provided as a Source Data file. Created in BioRender. Zhan, X. (2026) https:// BioRender.com/9bbcvwx.
In canonical GCs, B/plasma cells undergo SHM in dark zones (DZs) to generate mutation pools for clone selection and then migrate to light zones (LZs) for CSR alongside T-cell priming63. Since IgM+ PCZ and IgG+ PCZ exhibited exclusive expression of IGHM and IGHG genes (Fig. 4f, Supplementary Fig. 7b), we thus speculated that IgM+ PCZ and IgG+ PCZ might serve as the ectopic DZ-like and LZ-like structures with distinct clonal activities. First, we checked MKI67 (proliferation marker), ITGAX (dendritic cells, DC), LAMP3 (conventional DC3, cDC3)64,65, and CLEAC10A (conventional DC2, cDC2) expression of iTLSs and PCZs (Supplementary Fig. 8c, d). As expected, the expression of LAMP3 was highest in iTLSs41,61. IgM+ PCZs and IgG+ PCZs showed the lowest and highest expression of ITGAX and CLEC10A, respectively. The expression of MKI67 increased from inner iTLSs to IgG+ PCZs. The expression pattern of these genes by IgM+ PCZs and IgG+ PCZs phenotypically resembled the DZs and LZs in canonical GCs66,67,68. Next, we checked the IgH clone sharing of these two niches. Of note, IgM+ PCZ and IgG+ PCZ exhibited distinct IgH clone type numbers in different tumors (Supplementary Fig. 8e). LUAD_P1/P6/P11 had larger numbers of IgG+ PCZ IgH clone types, while other tumors showed larger numbers of IgM+ PCZ IgH clone types. We thus calculated the shared proportion using shared clone type numbers divided by IgM+ or IgG+ (the smaller one) PCZ-associated clone type numbers. As a result, IgM+ PCZ and IgG+ PCZ shared large proportions of clones (50 ~ 100%, Supplementary Fig. 8e), implicating that IgH clones in IgM+ PCZ and IgG+ PCZ might originate from the same B/plasma cell lineage precursors. By checking the expression of enzymes involved in SHM and CSR, we observed that IgM+ PCZ expressed the highest level of activation-induced cytidine deaminase (AICDA), an enzyme mediating SHM in DZs69,70, in line with the highest mutation frequencies in IgM+ PCZs (Fig. 5j). IgG+ PCZ expressed the highest level of uracil-DNA glycosylase (UNG) and apurinic/apyrimidinic endodeoxyribonuclease 1 (APEX1), enzymes mediating the CSR process in LZs71,72, as well as B-cell lymphoma 6 (BCL6) that marks class-switched B/plasma cells73(Fig. 5j). Furthermore, we examined the presence of T cells in PCZ. By quantifying T cells using the spatial transcriptome data, we observed a significantly higher proportion of T cells in IgG+ PCZ compared IgM+ PCZ (Fig. 5k). This result was further validated by the co-staining of CD3E/CD138. An abundant infiltration of T cells was observed in IgG+ PCZ but not in IgM+PCZ (Fig. 5l, m). The contrast T-cell proportions implicated a role played by T cells in IgG+PCZ, but not in IgM+PCZ. To explore this, we partitioned the PCZs using DBSCAN and obtained 113 spatially discrete PCZ clusters (Supplementary Fig. 8f). T cell signature genes in each PCZ cluster correlated with increased B/plasma cells pseudotime scores (p = 0.07, r = 0.516, Fig. 5n), CSR events (p = 0.04, r = 0.567, Fig. 5o), and IgG secretion (p = 0.005, r = 0.722, Fig. 5p), but not with the mutation accumulation of IgH chains (p = 0.17, Fig. 5q), suggesting a CSR-promoting role played by T cell in IgG+ PCZs but not in IgM+ PCZs. T cells are reported to participate in B-cell fate determination via secreting interleukins or through membrane ligand-receptor interaction74,75,76. We then plotted the potential molecular candidates expressed by T cells to mediate B/plasma cell effector function and differentiation. As shown by the matrix plot, IL6 was expressed in inner and outer iTLSs. Meanwhile, CD40LG was prominent in iTLSs and IgG+ PCZs (Supplementary Fig. 8g). These data imply that IL6 signals might be involved at a very early stage of intratumoral B-cell fate transition in iTLS. At the same time, sustained stimulation via CD40LG might drive the terminal differentiation of plasma cells in IgG+ PCZs.
Recent studies demonstrated the antigen-targeting cytotoxicity of germinal center-homing CD8+ T cells in certain infections77,78. We next investigated the tumor reactivity of the T-cell clones occurring in the IgG+ PCZ. According to the presence or absence of IgG+ PCZs, we stratified the TCRβ clones into present clones and absent clones. Notably, TCRβ clones present in IgG+ PCZs exhibited a pronounced aggregating pattern compared to those absent from these regions (Fig. 5r). To quantify this spatial clustering, we used Moran’s Index, a statistical measure of spatial autocorrelation based on Cartesian coordinate systems79. TCRβ clones located in the IgG+ PCZs showed significantly higher Moran’s Index values than IgG+ PCZ-absent clones, indicating enhanced spatial aggregation (Supplementary Fig. 8h). Antigen-driven clone expansion generates large populations of lymphocyte clones to mount antigen-targeting immune response80,81,82. Clone sizes, therefore, serve as a proxy for tumor reactivity. Compared to clones absent from IgG+ PCZ, those present in IgG+ PCZ consistently exhibited larger clone sizes across all 11 tumors (Fig. 5s), indicating robust clonal expansion. Furthermore, we analyzed the expression of effector function genes of these clones in the tumor region. In contrast, TILs of IgG+ PCZ-present TCRβ clones exhibited stronger killing effect and stemness represented by the expressions of tumor necrosis factor (TNF), T-cell factor 7 (TCF7), and Fas ligand (FASLG), and exhaustion markers (HAVCR2/PDCD1/TIGIT/LAG3) (Supplementary Fig. 8i). To assess the tumor reactivity of IgG+ PCZ-present and -absent clones, we applied two tumor reactivity-associated gene signature sets for scoring83,84. Clones present in the IgG+ PCZ consistently exhibited higher tumor-reactivity scores compared to their absent counterparts (Fig. 5t). Taken together, our data suggested that IgM+ PCZ are characterized by SHM, whereas IgG+ PCZ exhibit features of CSR and tumor-reactive T cell priming. These two spatial structures thus serve as ectopic DZ-like and LZ-like structures, respectively, collectively recapitulating the canonical GC reaction in the tumor (Fig. 5u).
Stereo-XCR-seq profiles exclusion of ectopic GC-like niches by terminally differentiated cancer-associated fibroblasts
Given the role of ectopic GC-like niches in priming tumor-reactive lymphocytes, we next investigated the factors potentially associated with B/plasma clone sizes and the clinical impact of IgG+ PCZs. Firstly, we quantified the top IgH clone sizes within each bin50 spot (Fig. 6a). As expected, >95% of non-expanded regions resided in tumors, whereas iTLS and PCZ harbored markedly larger clones (Fig. 6b). Interestingly, we observed that the IGHM expression concentrated in the regions with a clone size of two and three, while the expression of IGHG1/3/4 concentrated in the regions with a clone size over four (Fig. 6c). Consistently, bin50 spots with larger clone sizes exhibited co-localization with IgG+ PCZ, implicating an association between intensive B/plasma cells clonal expansion and IgG+ PCZ formation (Fig. 6d, e). Furthermore, we integrated PCZ, iTLS, myeloid cell region, and immune-excluded stroma into a single “stromal” meta-cluster and categorized these stromal spots into non-expanded regions (clone size <2) and expanded regions (clone size ≥2, Fig. 6a). The distinct spatial distributions of non-expanded and expanded stromal regions in the tumor (Fig. 6f) suggested that these two microenvironments possess differing characteristics. To explore which cells might be possibly involved in the B/plasma cell fate transitions, we further subdivided the major cell types into 40 clusters (Supplementary Fig. 9a). In addition to B cells, endothelial cells and plasma cells, we obtained 4 epithelial cell subclusters, 12 fibroblast cell subclusters, 7 myeloid cell subclusters and 14 T cell subclusters (Supplementary Fig. 9a and Fig. 6g). Pairwise cell-cell co-localization analysis showed that B cells co-localized with most T-cell subclusters (Fig. 6g). In contrast, plasma cells co-localized with CXCL13+ CD8 T cells, macrophages, and INHBA+ fibroblasts (Fig. 6g), which was further enhanced by the observations in all 11 patients (Fig. 6h). In addition, these three cell types exhibited increasing proportions in iTLS, PCZ, immune-excluded stroma and myeloid cell region (Supplementary Fig. 9b). In comparison to IgM+ PCZ, IgG+ PCZ exhibited higher proportions of these three cell types, indicating the associations towards the clone expansion of B/plasma cell (Supplementary Fig. 9c). Coinciding with the enrichment of tumor-reactive T cells in the ectopic LZ-like structures (Fig. 5k–t), we observed an enrichment of CXCL13+ CD8 T cells in the expanded stroma (Fig. 6i, left). In contrast, INHBA+ fibroblast and macrophage co-localized in non-expanded stroma (Fig. 6i, right), which was further validated by CD86 (macrophage intuitive marker) and PDGFRA (fibroblast intuitive marker) co-staining (Supplementary Fig. 9d). Consistent with previous studies showing that stromal cells and macrophages can orchestrate immune-suppressive milieus85,86,87,88,89, our observation indicated a potential linkage between the co-occurrence of macrophage-INHBA+ fibroblast and the exclusion of B/plasma cell clonal expansion.
Fig. 6: Stereo-XCR-seq profiles exclusion of ectopic GC-like niches by terminally differentiated cancer-associated fibroblasts.
a Graph shows the definition of non-expanded stroma and expanded stroma. Created in BioRender. Zhan, X. (2026) https:// BioRender.com/9bbcvwx. b Representative stacked bar chart of LUAD-P1 shows the frequencies of spatial spots with indicated clone sizes in different clusters, n = 11 samples. Bars are colored by spatial clusters as indicated above. c Dot plot shows the expression of the gene of interest in stroma spots at different clone sizes. Dots are colored by mean expression and sized by expression fraction. d Representative spatial plot of LUAD-P1 shows clone sizes of each stromal spots, n = 4896 cells, n = 11 samples analyzed, scale bar = 1 mm. Each dot represents a bin50 spot, colored and sized by clone sizes. e Representative spatial KDE plot of LUAD-P1 (n = 12,131 cells) shows the distributive density of IgG+ PCZ, n = 11 samples, scale bar = 1 mm. f Representative spatial KDE plot of LUAD-P1 shows the distributive density of non-expanded stroma (left, n = 4323 cells) and expanded stroma (right, n = 573 cells), n = 11 sample analyzed, scale bar = 1 mm. g Heatmap (left below) and dot plot (upper right) show the cell co-occurrence and marker gene expression. Dots are colored by mean expression and sized by expression fraction, n = 11 samples. h Heatmap showing the co-occurrence of plasma cells with other cell types across n = 11 samples. i Spatial plot shows the distribution of the cells of non-expanded stroma (left, n = 34,925 cells) and expanded stroma (right, n = 8946 cells), scale bar = 1 mm. Each polygon represents a cell, colored by cell subclusters. Representative FOVs in expanded stroma and non-expanded stroma are selected to present the co-localization of different cell clusters. N analyzed samples = 3 (LUAD-P1/P4/P11, with over 0.1% expanding stroma spots). j Dot plot showing the log10 fold change of the genes upregulated in expanded stroma spots. N analyzed samples = 3 (LUAD-P1/P4/P11, with over 0.1% expanding stroma spots). Genes are ranked by log10 fold changes. Each dot represents a gene; two-tailed unpaired Student’s t tests were used. No adjustments for multiple comparisons were applied. Genes of interest are listed above the plot, ****p < 0.0001. k Heatmap shows the expression of the gene of interest in stroma spots at different clone sizes. N analyzed samples = 3 (LUAD-P1/P4/P11, with over 0.1% expanding stroma spots). Each line represents a bin50 spot, ranked by clone size and colored by mean expression. l Survival test of LUAD patients (n = 478 samples) using IGHG1/3/4 and eight non-expanded stroma-related genes as a signature. Kaplan–Meier survival test, hazard ratio (HR) and significance tests (p value, p(HR)-values) are labeled. m Dot plot shows the expression of the gene of interest by fibroblasts, B/plasma cells and macrophages in scRNA-seq dataset (GSE207422)102, n = 15 samples. Dots are colored by mean expression and sized by expression fraction. Source data are provided as a Source Data file.
Next, we compared the gene expression profile between non-expanded stroma and expanded stroma to dissect the stromal programs opposing ectopic GC-like niche formation. Genes upregulated in the non-expanded stromal regions associated with fibroblast activation (CCL18)90, collagen production (COL10A1, COL5A1, COL12A1)91,92,93, ECM remodeling (MMP14, MMP11)94,95, and suppressive microenvironment establishment (IHNBA, TGFBI)96,97 (Fig. 6j). In addition to the co-expressing pattern of these genes in the non-expanded region (Fig. 6f, i), these genes exhibited a declining trend alongside the increased IgH clone sizes (Fig. 6k), revealing mutually exclusive localization of IgG+ PCZs and ECM-remodeling stroma. We then sought to identify the primary cellular sources of these genes. scRNA-seq data attributed these non-expanded stromal genes to INHBA+ fibroblasts and macrophage (Supplementary Fig. 10a), suggesting an ECM-promoting role of these two cell types in shaping the IgG+ PCZ-excluding microenvironment. In addition to the higher production of TGFBI and CCL18 by the macrophage in the non-expanded stroma, fibroblasts in the non-expanded stroma also exhibited a higher expression of ECM-related genes over their counterparts in expanded stroma or other spatial clusters (Supplementary Fig. 10b, c). Because macrophages could activate fibroblasts through CCL1890, we constructed a diffusion pseudotime to infer the fibroblast fate transition in the LUAD context98. By setting starting point at CD34+ fibroblast, which exhibited distinct stemness among all fibroblast subsets99 (Supplementary Fig. 10d), we observed two terminally differentiated status of the intratumoral fibroblasts (Supplementary Fig. 10e). In addition, two potential fate transition paths constructed by PAGA pointed at the expanded regions (path1) and the non-expanded regions (path2), respectively (Supplementary Fig. 10f), coinciding with the advanced pseudotime of the fibroblasts in the non-expanded regions (Supplementary Fig. 10g). Given the terminally differentiated fibroblasts from the non-expanded regions upregulated genes related to ECM remodeling (Fig. 6j, k), we ran gene ontology (GO) term analysis using genes significantly upregulated (p < 0.05) in non-expanded regions and expanded regions, respectively. In contrast to the upregulated BCR signaling in expanded regions, the non-expanded regions exhibited enhanced activities of collagen fibril organization and extracellular structure organization (Supplementary Fig. 10h). Together, these data implicate that TGFBI+macrophages might promote ECM organization via CCL18-mediated promoting fibroblast activation to exclude lymphocyte clone expansion.
Since IgG+ PCZ exhibited tumor-reactive clone selection and expansion, we investigated the clinical impact of IgG+ PCZ formation. We calculated IgG+ PCZ signature scores (Z-scoring IGHG1/3/4) and IgG+ PCZ-excluding signature scores (Z-scoring genes upregulated in the non-expanded stroma, highlighted in Fig. 6h) using the bulk RNA-seq expression matrix of the 478 LUAD tumors from the TCGA database100,101. By integrating the survival information, we observed that patients with higher IgG+ PCZ scores exhibited improved overall survival (Fig. 6l, left), whereas patients with higher IgG+ PCZ-excluding scores exhibited poorer overall survival (Fig. 6m, right), highlighting the prognostic relevance of IgG+ PCZ formation in LUAD patients. Furthermore, we examined gene expression profiles of 1122 fibroblasts, 18,501 macrophages, and 3661 B/plasma cells from an ICB-treated LUAD scRNA-seq dataset (GSE207422, 15 patients)102. Expression of MMP11, MMP14 and COL12A1 by fibroblast and TGFBI by macrophage in pre-treatment tumors associated with marginal disease regression (stable disease, SD), whereas expression of IGHG1/2/3/4 by B/plasma cells associated with improved response (partial response, PR, Fig. 6m). In an additional scRNA-seq dataset (GSE243013) from post-ICB treatment LUAD103, we also observed that ICB responders exhibited the highest clonal expansion of CXCL13+ CD8 T cells (Supplementary Fig. 10i, upper panel), whereas non‑responders were characterized by the highest clonal expansion of ANXA1+ CD4 memory T cells (Supplementary Fig. 10i, lower panel). This pattern well aligns with our observation that CXCL13+ CD8 T cells represented a major tumor‑reactive T‑cell subtype that co-localizes with IgG+ plasma cells in ectopic LZ-like structures. Collectively, our data suggested that the pre-existence of ECM-modulating macrophage and fibroblasts in LUAD may be prognostically unfavorable and correlating with reduced immunotherapy responsiveness and overall survival. In contrast, the formation of IgG+ PCZ is associated with the expansion of tumor-reactive clones, indicating a potentially improved response to ICB treatment.
