A single-cell atlas of non-haematopoietic cells in human lymph nodes and lymphoma reveals a landscape of stromal remodelling

The activities of non-haematopoietic cells (NHCs), including mesenchymal stromal cells and endothelial cells, in lymphomas are reported to underlie lymphomagenesis. However, our understanding of lymphoma NHCs has been hampered by unexplained NHC heterogeneity, even in normal human lymph nodes (LNs). Here we constructed a single-cell transcriptome atlas of more than 100,000 NHCs collected from 27 human samples, including LNs and various nodal lymphomas, and it revealed 30 distinct subclusters, including some that were previously unrecognized. Notably, this atlas was useful for comparative analyses with lymphoma NHCs, which revealed an unanticipated landscape of subcluster-specific changes in gene expression and interaction with malignant cells in follicular lymphoma NHCs. This facilitates our understanding of stromal remodelling in lymphoma and highlights potential clinical biomarkers. Our study largely updates NHC taxonomy in human LNs and analysis of disease status, and provides a rich resource and deeper insights into LN and lymphoma biology to advance lymphoma management and therapy. Abe et al. profile, characterize and compare non-haematopoietic cells in normal human lymph nodes versus nodal lymphomas from patients, providing insights into stromal modelling in health and disease.

L ymphomas are haematological malignancies that often develop from LNs 1 . Despite advances in treatments, most lymphoma subtypes remain incurable. Therefore, new therapeutic approaches are needed, including those that target the tumour microenvironment 2,3 . In lymphomas, as in solid cancers 4,5 , the activities of NHCs, such as mesenchymal stromal cells (SCs) and endothelial cells, are thought to facilitate lymphomagenesis and therefore have potential as therapeutic targets 2,3 . Indeed, some lymphoma subtypes are reported to exhibit unique interactions with NHCs [6][7][8][9] ; however, lymphoma NHC research is far behind that of solid cancers 10 . In particular, follicular lymphoma (FL) cells are considered to actively interact with NHCs to promote emergence and expansion 9,11,12 . SC-derived CXCL12 recruits FL cells in cooperation with CXCL13 (produced by follicular dendritic cells (FDCs)), which contribute to the follicular localization of tumour cells and their proliferation 9,13 . Other FDC-derived molecules, including BAFF (encoded by TNFSF13B), interleukin-15 and HGF, may have anti-apoptotic effects on FL cells [14][15][16] . Unfortunately, a complete understanding of the temporal and spatial associations that underlie these activities is hampered by the heterogeneity of NHCs. In fact, definitive NHC classification has not yet been achieved in humans, even in normal LNs 17,18 . Moreover, the identification of alterations in LN NHC (LNNHC) heterogeneity in the context of lymphomas is barely underway.
To address this issue, scRNA-seq was used in this study to construct an atlas of human NHCs in LNs and lymphoma. We aimed to identify previously unrecognized NHC heterogeneity in human LNs and to distinguish NHCs from lymphomas to define the global influences of lymphoma cells on the NHC niche. This approach can provide deep insights into lymphoma stromal biology and resources applicable to future studies of lymphomas and identify potential stroma-derived biomarkers that may serve as clinical indicators and/or therapeutic targets. Table 5), which is upregulated following induction of cell-cycle progression 33 , and SELE (Fig. 2h,i and Extended Data Fig. 2a,b), which is upregulated by inflammation 34 . By contrast, hHEVs expressed SELE at low levels ( Fig. 2i and Extended Data Fig. 2a,b). Notably, C-aHEVs and aHEVs both expressed stress-related genes, including those associated with heat shock proteins, NF-κB activation, JNK activation and shear stress 35 (Extended Data Fig. 2b  and Supplementary Table 6), which suggests that these subclusters respond to active cell deformation or damage.
We next performed trajectory analysis on integrated MFLN BEC data using the Monocle 3 pipeline 36 . We were able to identify all BEC subclusters in a cell object generated using Monocle 3 (Extended Data Fig. 2c). Trajectory of the arterial component flowed from ABECs to aBECs and cBECs, finally reaching tBECs (Fig. 2j). Similarly, trajectory of the venous component initially traced HEV subclusters (aHEVs and hHEVs), then proceeded to capillary subclusters (C-aHEVs and cBECs) and finally to tBECs (Fig. 2j). These findings support the idea that tBEC migration in LNs generates new capillary BECs 37 .
Gene ontology (GO) analysis revealed that factors involved in blood vessel development are enriched in ABECs, caBECs, aBECs and tBECs ( Fig. 2k and Supplementary Table 7), which is in agreement with their arterial or tip cell characteristic. Leukocyte migration and cellular extravasation signatures were most enriched in aHEVs (Fig. 2k). Molecules associated with apoptosis were enriched in C-aHEVs and aHEVs (Fig. 2k). Moreover, as reported in mice 25 , CXCL10-HEVs expressed molecules associated with interferon and cytokine signalling (Fig. 2k).
Immunofluorescence (IF) analysis of BECs stained with GJA5, SSUH2 or INSR identified them as large arterial BECs in LNs (ABECs), arterial BECs outside LNs (caBECs) and arterioles (aBECs), respectively (Extended Data Fig. 2d-f). We also detected tBECs as cells stained positive for LY6H or PGF in the tips of PLVAP + cBECs ( Fig. 2l and Extended Data Fig. 2g). HEVs strongly expressing SELE (aHEVs) were frequently observed in interfollicular regions (IFRs) (Fig. 2l), which indicates that IFRs may serve as niches that play pivotal roles in promoting the influx of immune cells into LNs. Moreover, staining for PLVAP, HES1 and the HEV marker MECA-79 revealed that PLVAP + HES1 + capillaries (C-aHEVs) and MECA-79 + HES1 + HEVs (aHEVs) (Extended Data Fig. 2b) were localized near each other (Extended Data Fig. 2h). Notably, CXCL10-HEVs were frequently observed in IFRs and were localized exclusively in the vicinity of aHEVs ( Fig. 2l and Extended Data Fig. 2i). These findings, together with the GO analysis, suggest that rare CXCL10-HEVs may activate cellular trafficking of adjacent HEVs through cytokine signalling, which results in the heterogeneity of human HEVs.
In summary, our single-cell atlas of LN BECs identified three, three and four transcriptionally distinct subclusters in arterial, capillary and venous BECs, respectively, which demonstrates the unique heterogeneity of these cells in humans (Extended Data Fig. 2j-l).
Twelve subclusters of human LN NESCs. NESCs were divided into the following 12 subclusters: SCs at capsule adventitia (advSCs); SFRP4 + SCs (SFRP4-SCs); SFRP2 + SCs (SFRP2-SCs); SCs enriched for tumour necrosis factor (TNF) signalling (TNF-SCs); C7 + SCs (C7-SCs); AGT + SCs (AGT-SCs); TRCs; pericytes (PCs); smooth muscle cells (SMCs) with high or low ATF3 expression (ATF3 hi or ATF3 lo SMCs); MRCs; and FDCs ( Fig. 4a,b and Extended Data Fig. 4a). TRCs, PCs, MRCs and FDCs were annotated on the basis of conventional taxonomy 18 . The results that contributed to the annotation of these known subclusters are included in the Supplementary Note. DEG analysis revealed that NESC subclusters exhibited more than 100 DEGs each (Fig. 4c). advSCs showed the highest CD34 expression among NESCs (Fig. 4d) and are considered the human counterpart of mouse CD34 + SCs observed at adventitia of the LN capsule 21 . Both SFRP4-SCs and SFRP2-SCs shared SFRP2 expression and were discriminated by higher SFRP4 expression in the former (Fig. 4d,e). SFRP4-SCs also showed relatively high INMT expression (Supplementary Table 10), which suggests that they are the counterpart of mouse Inmt + SCs observed exclusively at medullary cords 21 . TNF-SCs were specifically characterized by PTX3 expression, and C7-SCs by abundant C7 expression (Fig. 4d,e). AGT-SCs expressed AGT and high levels of the apolipoprotein genes APOE and APOC1 (Fig. 4d,e). ATF3 hi and ATF3 lo SMCs both expressed muscle-specific MYH11 and PLN (Fig. 4d,e), but differed in the expression of genes associated with cellular responses to stress or mechanical stimuli (Fig. 4f and Supplementary Table 11). Notably, TNFSF13B, which encodes B-cell-activating factor belonging to the TNF family (BAFF) and is thought to define FDCs 20 , was expressed by both MRCs and FDCs but at higher levels by MRCs (Fig. 4d,e).
Trajectory analysis revealed that MRCs were connected to TNF-SCs and C7-SCs ( Fig. 4g and Extended Data Fig. 4b), which indicates that the latter two subclusters might differentiate into MRCs. Additional analysis showed a continuous trajectory from SMC subclusters to PCs, TRCs, MRCs and finally to FDCs in human LNs ( Fig. 4h and Extended Data Fig. 4c), which is consistent with findings in mice of fibroblastic reticular cells in the splenic white pulp 38 .
GO analysis revealed that advSCs expressed high levels of genes involved in the formation of elastic fibres and the extracellular matrix (ECM) 39 (Fig. 4i and Supplementary Table 12). In agreement with the preferential localization of mouse Inmt + SCs at the medulla 21 , their human counterparts, SFRP4-SCs, expressed high levels of genes involved in ECM formation ( Fig. 4i and Supplementary  Table 12). TNF-SCs expressed genes associated with TNF signalling (IL6 and CCL2) ( Fig. 4i and Supplementary Tables 10 and 12), which suggests that they function in the chemotaxis of CCR2-expressing T cells, monocytes and dendritic cells to antigen sites 40 . C7-SCs expressed genes related to chemotaxis regulation (Fig. 4i), including CXCL12 (Fig. 4d), which supports transendothelial T-cell migration across HEVs 41 . Top DEGs for AGT-SCs included APOE, AGT and LPL, which participate in remodelling of protein-lipid complexes and plasma lipoprotein particles ( Fig. 4i and Supplementary Tables 10 and 12), which suggests that AGT-SCs may participate in lipid metabolism or transport.
IF staining was performed to identify the localization of each subcluster in LNs (Extended Data Fig. 4d-o). Fibroblasts positive for decorin (encoded by DCN), a strong marker of advSCs, SFRP4-SCs, SFRP2-SCs, TNF-SCs and C7-SCs (Extended Data Fig.  4d), were widely distributed in the adventitia, IFRs and medulla (Extended Data Fig. 4e). FBN1 + SCs (advSCs) (Fig. 4d) were observed at the capsule adventitia, as observed in mice 21 (Extended Data Fig. 4f). SFRP2 + SCs (SFRP2-SCs and SFRP4-SCs) were preferentially distributed in the medulla (Extended Data Fig. 4g). PTX3 + SCs (TNF-SCs) were observed in IFRs (Extended Data Fig. 4h). C7-SCs were most frequent in the outer cortex, excluding follicles (Extended Data Fig. 4i), which is consistent with their proposed role in facilitating immune cell migration. AGT + cells were found on outer regions of the IFRs, frequently situated between SCSs and HEVs (Extended Data Fig. 4j). SMCs were observed as α-smooth muscle actin + (encoded by ACTA2), MYH11 + or PLN + cells (Fig. 4d) around not only arterial BECs but also some HEVs (Extended Data Fig. 4k,l). ATF3 was positive in some SMCs around HEVs in the IFRs (aHEVs), as well as around arteries (Extended Data Fig. 4l). In line with the DEG analysis between SMC subclusters, ATF3 + SMCs were also marked by HSP70 (encoded by HSPA1A) expression (Extended Data Fig. 4m), which probably reflects cell damage induced by blood flow 42 and/or immune cell trafficking.
To summarize, we identified 12 NESC subclusters, thereby showing unanticipated heterogeneity, linked to the distribution of other NHC subsets and LN niches (Extended Data Fig. 4p-r).
We accomplished a single-cell atlas of NHC components in human LNs (Extended Data Fig. 4p). Additional   Remodelling of NHC proportions in FL. Using this atlas, we next sought to explore alterations in FL NHCs at subcluster levels by comparing them with MFLN counterparts (Fig. 5a and Extended Data Fig. 7a). Overall, the proportion of BECs were markedly increased in FL relative to MFLNs, whereas the proportion of LECs decreased (Fig. 5a). Moreover, the proportion of arterial subclusters were increased in FL BECs (Fig. 5a). In FL NESCs, the proportion of FDCs was greatly increased (Fig. 5a). Notably, MRCs were also greatly increased in FL, whereas advSCs, SFRP4-SCs, SFRP2-SCs and TNF-SCs were decreased (Fig. 5a).

ResouRce
NaTurE CELL BiOLOgy some FL NESC subclusters, although this finding was not significant (Extended Data Fig. 7b). Notably, in some NESC subclusters, we observed marked upregulation of genes relevant to solid cancers but previously not associated with lymphomagenesis. Among them, POSTN, which encodes periostin (a protein secreted by cancer-associated fibroblasts (CAFs) and promotes the formation of cancer stem cell, perivascular and premetastatic niches 43 ), was substantially upregulated in TRCs and PCs and in SMC subclusters of FL (Fig. 5c, Extended Data Fig. 7b and Supplementary Table 17). The expression of EGFL6, which encodes EGFL6 (a member of the EGF-like superfamily that reportedly promotes tumour cell growth by stimulating angiogenesis 44,45 ), was highly upregulated in TRCs, SMC subclusters and MRCs (Fig. 5c 49,50 . Indeed, FL HEV subclusters expressed genes that regulate cellular adhesion and migration (Extended Data Fig. 7c and Supplementary Tables  15 and 18). Notably, expression of the tip cell markers LY6H, PXDN, PGF and LOX was markedly upregulated in FL tBECs ( Fig. 5c and Supplementary Table 15), which suggests that they are involved in the acceleration of angiogenesis. The significant decrease in the proportion of LECs in FL suggests that there is widespread lymphatic damage. IF staining confirmed that the LEC density was lower in FL compared with that in MFLNs (Extended Data Fig. 7d,e). Many FL LEC subclusters also showed upregulation of heat shock genes as well as CD74, which reportedly functions in wound healing 51 (Fig. 5c, Extended Data Fig. 7c and Supplementary Tables 16 and  18). CD74 overexpression was confirmed in FL LECs by IF staining (Extended Data Fig. 7f,g).

Landscape of intercellular interactions in FL stroma.
To assess the NHC-malignant B-cell crosstalk underlying FL growth, we performed scRNA-seq of cryopreserved CD45 + cells from nine FL samples (FL 2-FL 10) and extracted gene expression profiles of malignant B-cell clusters in silico from each (Extended Data Figs. 1c and 8a-d). We then performed intercellular ligand-receptor interaction analyses between FL NHC subclusters and malignant B cells using CellPhoneDB 52 . Thereafter, we extracted significant interactions that were considered upregulated in FL NHC subclusters relative to the corresponding MFLN subclusters.
We identified a total of 58 interactions, including some previously uncharacterized in FL (Fig. 6a). In BECs, we noted that overexpression of JAG1, which is reportedly observed in B-cell lymphoma BECs and associated with aggressive lymphoma phenotypes 53 , was limited to only larger arterial BEC subclusters (ABECs and caBECs) (Fig. 6a). Interactions mediated through adhesion molecules, including the SELE-CD44 interaction 54,55 , were activated mainly in HEV subclusters (C-aHEVs, aHEVs and hHEVs) (Fig. 6a), which suggests that these HEV subclusters may contribute to the haematogenous expansion of FL cells 54,56 . Interactions that promote cancer cell death and mediated by TNFSF10 were markedly upregulated in several LEC subclusters 57 (Fig. 6a), which suggests that LECs may antagonize lymphoma development. In NESCs, interactions associated with TNF signalling, cell adhesion, PDGF signalling and chemokine signalling were differentially activated among subclusters (Fig. 6a). Notably, overexpression of CXCL12, which reportedly supports FL cell migration, adhesion and activation 58 , was observed in advSCs (Fig. 6a). Moreover, interactions via BAFF were upregulated, even in medullary SCs (SFRP4-SCs), which suggests that stromal remodelling in FL supports the extrafollicular expansion of malignant B cells 59 . In advSCs and medullary SC subclusters, interactions mediated by stroma-derived CD70 were enhanced (Fig. 6a). Interactions mediated through PDGFRB, which promotes cell migration and angiogenesis 60 , were enhanced in TRCs and PCs (Fig. 6a), which suggests that during FL expansion, mechanisms other than CCR7-CCL19/CCL21 signalling may drive the homing of malignant B cells to the T-cell zone 61 . Instead, the CCR7-CCL19 interaction was extended to non-TRC SCs (TNF-SCs and PCs) (Fig. 6a). Consistent with the DEG analyses of MFLNs and FL, the CXCL13-CXCR5 axis 9,13 was activated in MRCs and FDCs (Fig. 6a).
Enhanced CD70-CD27 interaction across FL stroma. Based on the above interactome analysis results, we next sought to explore an interaction that can potentially be targeted in lymphoma. We carefully surveyed candidate interactions from the perspective of novelty in the field. We noted that the CD70-CD27 interaction in solid and haematological cancers has attracted increasing attention [62][63][64][65] , whereas interactions mediated by stroma-derived CD70 have rarely been investigated. Accordingly, we focused on the CD70-CD27 interaction for functional validation to verify the usefulness of our atlas-based analyses and to propose a potential mechanism in the stroma relevant to FL progression. Initially, we confirmed that

NaTurE CELL BiOLOgy
CD70 is overexpressed in FL medullary and adventitial SCs by IF staining (Fig. 6b,c). We next examined the gene and protein expression levels of the CD70 ligand CD27 in the B cells of FL samples. Single-cell transcriptomic analysis of FL B cells showed that CD27 was significantly upregulated in malignant B cells compared with non-malignant B cells (Extended Data Fig. 8e). Consistent with these results, flow cytometry analysis of FL haematopoietic cells showed that the CD19 + CD10 + cell population (malignant B-cell enriched fraction) in 5 out of 8 (62.5%) biologically independent samples was positive for CD27, and its expression was also significantly higher in the CD19 + CD10 + population than in the CD19 + CD10 − population (non-malignant B-cell fraction) (Extended Data Fig. 8f,g). Among the five CD27 + FL samples, four (80.0%) showed unequivocal binding to recombinant human CD70-Fc protein (Fig. 6d). The binding of malignant B-enriched cells to CD70-Fc protein was significantly inhibited by the treatment of the cells with an anti-CD27 function-blocking antibody in all four cases (Fig. 6d,e). Next, we performed ex vivo cell adhesion assays using    ResouRce NaTurE CELL BiOLOgy significantly decreased following treatment with the anti-CD27 antibody (Fig. 6f,g).

Prognostic implications of stroma-derived markers in FL.
Next, we tested the applicability of our single-cell analysis of NHCs in the search for clinically relevant factors. To correlate niche-specific or subcluster-specific alterations in NHCs with survival of patients with FL, we utilized a bulk microarray dataset of 180 FL biopsy samples from newly diagnosed patients with available survival information 66 .
We narrowed down multivariate analysis candidates to seven genes (LY6H, LOX, PTGIS, TDO2, REM1, PIEZO2 and CHI3L1) expressed at minimal levels in FL haematopoietic cells but at high levels in FL BEC or NESC subclusters compared with MFLN counterparts. We hypothesized that they were probably associated with unfavourable prognosis (Fig. 7a and Extended Data Fig. 9a-e).
In the multivariate analysis, increased expression of the tip cell markers LY6H and LOX, as well as TDO2 and REM1, were associated with an unfavourable prognosis, even after adjustment for the international prognostic index 67 (Fig. 7b).
For each of the four genes, we performed IF staining in MFLN and FL samples. Cells expressing LY6H, LOX, TDO2 or REM1 were increased in FL compared with those in MFLNs (Fig. 7c).
Findings from the additional prognostic analyses are described in the Supplementary Note, Extended Data Fig. 9f,g and Supplementary Table 19.
Observation of NHC subclusters across lymphomas. Finally, we examined whether our single-cell atlas was applicable to different lymphoma subtypes. We also investigated a more aggressive FL stromal remodelling phenotype. To this end, we performed scRNA-seq of stroma-enriched cells from five nodal peripheral T-cell lymphoma (PTCL) samples and three diffuse large B-cell lymphoma transformed from FL (tDLBCL) samples (Fig. 1a, Extended Data Fig. 1c and Supplementary Table 20). Detailed findings are described in the Supplementary Note and Extended Data Fig. 10a-h. In brief, NHC subclusters were detectable in PTCL and tDLBCL data and we observed distinct alterations in tDLBCL stroma that probably represented a terminal form of stromal remodelling in FL.

Discussion
Here we presented a human LNNHC map at single-cell resolution that was useful for exploring changes in lymphoma NHCs.
First, we shed light on differences in mouse and human LNNHC heterogeneity. Overall, our findings suggest that human LNs harbour unique NHC subpopulations that have not been detected in murine LNs, which emphasizes the need for further human studies. As in previous studies that used fresh human LN samples 22,26 , we used MFLNs from patients with tumours to construct the atlas. Notably, a detailed analysis of LNNHCs from an individual with a benign tumour indicated that the clustering was comparable between samples from the individual with a benign tumour and patients with a malignant tumour (Supplementary Note). This observation suggests minimal or negligible influence of malignancy-derived factors on our atlas.
Second, multistep DEG analyses revealed subcluster-specific changes in FL, including those with a previously unknown function in lymphoma. We found that upregulation of some of the known intercellular interactions across FL NHCs and malignant B cells extended to unanticipated NHC subclusters and, conversely, other interactions were enhanced in limited NHC subclusters. These observations largely increase the resolution of our understating of stromal remodelling in lymphoma. Additionally, these findings may be of clinical importance, as these may be considered as potential stroma-derived prognostic factors. Notably, two tip-cell markers were upregulated in FL and could serve as prognostic factors. LOX enzymatic activity is reported to drive tumour angiogenesis by activating PDGFRβ signalling in vascular SMCs, which is consistent with our DEG analysis that FL SMCs expressed PDGFRB at high levels 68 . Meanwhile, our observation of LY6H expression in tip cells has not previously been described in mouse or human endothelial cells. We also identified TDO2 as a prognostic predictor of FL. TDO2 may function to attract regulatory T cells, antagonize CD8 + T-cell activity and accelerate myeloid cell tolerogenicity 69 . REM1 overexpression in TRCs and PCs was also associated with unfavourable FL outcomes. Thus, further analysis of this gene, which has been scarcely explored, and its relevance to the lymphoma stroma is warranted. As the enrichment of the FL TRC signature per se was not prognostic (Supplementary Note), qualitative rather than quantitative alterations in certain NHC subpopulations may affect the chemoresistance and prognosis of FL more precisely. In addition to these prognostic factors, many upregulated genes with or without a known pro-tumorigenic function were included in our dataset, which makes our atlas a powerful discovery tool for additional therapeutic targets.
Third, we found that the CD70-CD27 interaction via stroma-derived CD70 was enhanced in FL. Although the role of CD70 has increasingly been investigated in the context of interplays across various immune cells and cancer cells 64,65 , lymphoma SCs have not been explored as a source of CD70. A recent report suggested that CD70 expressed by CAFs supports tumour progression in solid cancers by facilitating cancer cell migration 70 . Consistent with these findings, we confirmed binding between CD70 and malignant B cells that could be blocked by an antagonist against the CD70 ligand CD27. CD70 was upregulated in extrafollicular FL SCs, which suggests that CD70 may facilitate the infiltration of lymphoma cells into extrafollicular regions during tumour progression. Our analysis therefore proposes stroma-derived CD70 as a potential biomarker and therapeutic target for FL.
Last, we found that NHC heterogeneities in LNs were detectable even in aggressive lymphomas, thereby confirming the usefulness of our NHC atlas to characterize the stroma of various lymphoma subtypes. In particular, alterations in tDLBCL stroma harmonized with those in FL, thereby supporting the findings from the analysis of FL stroma. Furthermore, our findings indicate that extrafollicular SCs, including TRCs and medullary SCs, not only promote extrafollicular infiltration of FL cells but simultaneously differentiate into FSCs and are finally replaced by FSCs in more advanced phenotypes. This reflects a unique stromal transition that corresponds to the FSC-dependent growth of FL 71 .
Limitations of this study include the quantity of samples, which may not be sufficient to identify all NHC subpopulations or to precisely determine the correlation between the NHC heterogeneities in the transcriptome data and patient characteristics, such as genomic alterations. Second, we cannot completely exclude the possibility that our atlas is influenced by unknown factors from a distant malignancy. Third, our study was not designed to analyse other lymphoma subtypes or non-lymphoma diseases. Finally, further functional validation is required to confirm our findings relevant to each NHC subcluster.
In summary, our LNNHC atlas is of value to lymphoma researchers as it largely updates the NHC taxonomy in human LNs in the context of lymphoma. This study provides a platform for future research that aims to deepen our understanding of LN or lymphoma biology and to improve lymphoma management.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41556-022-00866-3.

Human samples. This study was approved by the Ethics Committee of the University of Tsukuba Hospital and the review boards of associated institutions that provided human samples (Kameda Medical Center, NTT Medical Center
Tokyo and Mito Medical Center) and conducted according to all relevant ethical regulations regarding human participants. Written informed consent was obtained from all participants. The participants were not compensated for taking part in the study. For scRNA-seq, MFLN samples were prospectively collected from patients with a neoplasm (n = 9) who had undergone surgical LN dissection between January and June 2020. Non-sentinel LNs without enlargement (<1 cm) were used. The collected LNs were verified as malignancy-free via flow cytometry analysis of pan-cytokeratin negativity (Extended Data Fig. 1b). Nodal FL (n = 10), PTCL (n = 5) and tDLBCL (n = 3) samples were also prospectively collected between August 2019 and May 2020. Furthermore, for functional experiments, additional nodal FL samples (n = 8) were collected between May 2020 and August 2021. Lymphoma diagnosis of tissue specimens was made pathologically, phenotypically and/or referring to results of cytogenetic examinations, including fluorescence in situ hybridization analysis by expert haematopathologists.
Library preparation, sequencing and data pre-processing. Sorted CD45cells were converted to barcoded scRNA-seq libraries using Chromium Single Cell 3′ reagent kits (V3) (10X Genomics) according to the manufacturer's instructions (CG000183 Rev A), aiming for 5,000-8,000 cells per library. Library quality control and quantification were performed using a KAPA Library Quantification kit for Illumina platforms (Kapa Biosystems, KK4873) and a 2100 Bioanalyzer High Sensitivity DNA kit (Agilent, 5067-4626). Libraries were sequenced on an Illumina HiSeq X Ten system with an average depth of 31,439 reads per cell, then mapped to the human genome (build GRCh38) and demultiplexed using CellRanger pipelines (10X Genomics, v.3.1.0).

Data processing and cell clustering of individual cases.
Pre-processed data from each sample were further processed and analysed individually using the R package Seurat (v.3.2.2) on RStudio (v.3.5.0 or v.4.0.2). After removing ribosomal genes, genes expressed in fewer than 3 cells and cells expressing fewer than 200 genes, we filtered out cells with fewer than 200 unique feature counts (low-quality cells). Cells with unique feature counts greater than three times the median value (possible doublets) and/or cells with more than twice the median number of mitochondrial genes (possible apoptotic or lysed cells) were also removed. We then normalized data using the NormalizeData function and extracted highly variable features using the FindVariableFeatures function. Normalized data underwent a linear transformation (scaling) and principal component analysis (PCA) based on variable features using the RunPCA function. Graph-based clustering was then performed according to gene expression profiles using the FindNeighbors and FindClusters functions with default parameters, and results were visualized using a nonlinear dimensional reduction UMAP technique running RunUMAP and DimPlot functions. Cell clusters were annotated based on the expression of canonical markers, including PECAM1 and JAM2 for BECs, PECAM1 and PROX1 for LECs, ACTA2 for SMCs, CCL19 and CCL21 for TRCs, CR2 for FDCs, DCN for other NESCs, PTPRC for contaminating lymphocytes, SDC1 for plasma cells, and CCR7 and CD83 (in cells weakly PTPRC + ) for dendritic cells 22,72,73 . MKI67 and TOP2A expression levels were used to identify clusters of an aggressively proliferative nature. For the FL sample 3, which came from a patient with intra-submandibular gland FL, any distinct clusters negative for all canonical markers and positive for keratin genes (indicating glandular tissue contamination) were removed. All other cases were confirmed to consist solely of these major clusters. We confirmed a negligible presence of ambient RNA contamination in single-cell NHC data, and found an imperceptible influence of potential RNA contamination on clustering results in all LN and lymphoma samples by using the DecontX (in the celda package, v.1.6.1) and SoupX (v.1.5.2) packages (data not shown) 74,75 .

Data integration with batch effect collection.
We performed canonical correlation analysis 76 to identify shared sources of variation across multiple datasets using the FindIntegrationAnchors function and integrated them using anchors from the IntegrateData function with canonical correlation dimensions of 20. Integrated data were scaled and underwent PCA as performed in individual datasets.
Supervised annotation and unsupervised clustering of LNNHCs. We performed graph-based clustering of PCA-reduced integrated data and supervised annotation, as described in 'Data processing and cell clustering of individual cases' above. Clusters characterized by extremely low unique feature counts (low-quality cells) were removed.
Next, we extracted the three major NHC components (BECs, LECs and NESCs) in silico and performed scaling, PCA-based dimensional reduction and unsupervised graph-based subclustering of each component. We removed subclusters that were considered possible doublets as characterized by high expressions of marker genes for different NHC components and incongruously high unique feature counts. In BEC subclustering, we also performed supervised annotation for the identification of arterial, capillary and venous BECs using canonical markers for each BEC component [29][30][31][32] . DEG analysis. DEG analysis was performed using the FindMarkers or FindAllMarkers functions with a minimum of 20% of the gene-expressing cells and a minimum log fold-change of 0.25 in gene expression between each cluster and other clusters. We primarily used the Wilcoxon rank-sum test for DEG detection. To confirm detected DEGs, we also used the model-based analysis of single-cell transcriptomics (MAST) method 77 . DEGs were defined as genes confirmed to show an adjusted P value (based on the Bonferroni correction) of <0.05 by using both methods. Results of the Wilcoxon rank-sum test were used to construct DEG lists and volcano plots. Volcano plots were created using the R package EnhancedVolcano (v.1.8.0). DEG analysis to compare corresponding clusters between mLN and pLN samples and between MFLN and FL samples was performed in a similar manner using the cut-off parameters described above.
For DEG analyses between MFLN and FL NHC subclusters, we adopted a multistep approach. Several previous studies had indicated differences in gene and protein expression between mLNs and pLNs [78][79][80][81] . Therefore, we initially profiled DEGs between mLNs and pLNs among MFLNs at subcluster levels (Supplementary  Table 13 and Supplementary Note). Referring to this profile, we identified DEGs upregulated in FL by removing those detected between mLNs and pLNs. We also performed DEG analysis between MFLN and FL NHC subclusters using only pLN samples (MFLN 7-MFLN 9 and FL 2-FL 10) to support the reliability of the detected DEGs.

Trajectory analysis.
We performed trajectory analysis using the Monocle 3 package (v.0.2.3) 36 in RStudio on integrated BEC, NESC and LEC data constructed using Seurat. Data pre-processing was performed using the preprocess_cds function, with the number of dimensions set at 100. Dimensionality reduction and clustering were performed using the reduce_dimension and cluster_cells functions, respectively. We then fit a principal graph within each cluster using the learn_graph function and visualized the order of cells in pseudotime by plot_cells or plot_cells_3d functions, as appropriate with the pseudotime colouring option.

Single-cell analysis of FL haematopoietic cells.
We performed single-cell analysis of cryopreserved CD45 + cells from nine FL samples (FL 2-FL 10). After thawing, cell suspensions were filtered through a 70-μm mesh and incubated with 7-AAD viability staining solution for 10 min in the dark. The 7-AADlive cells were sorted using a FACSAria II or III after removing doublets, then were converted to barcoded scRNA-seq libraries, as performed for CD45cells. Library preparation, sequencing and data processing were performed as for CD45cells. Data quality control, processing and graph-based clustering were performed in each individual case using the Seurat package, with dimension and resolution parameters of 50 and 0.5, respectively. Thereafter, we identified malignant B-cell populations by detecting restrictions of light chain kappa/lambda genes, as suggested by previous studies 83,84 . In brief, we projected the B-cell marker CD79A and the light chain genes IGKC (for light chain kappa) and IGLC2 (for light chain lambda) to cell clusters on the UMAP plot of each sample (Extended Data Fig. 8a). We then calculated the ratio of cells expressing IGLC2 and IGKC with expression levels of >1 and >2, respectively, in each B-cell cluster. We defined B-cell clusters with a ratio of >2.0 or <0.25 as malignant (Extended Data Fig. 8b).

NaTurE CELL BiOLOgy
Malignant B-cell signature analysis in FL B cells. To support the reliability of malignant B-cell detection, we performed signature analysis on data from FL B cells. We developed a gene set that represents a malignant B-cell signature based on the recent single-cell analysis of FL B cells reported by Andor et al. 83 . We carefully selected genes that were described as significantly upregulated in malignant compared to non-malignant B cells in a uniform manner among different FL samples 83 . Selected genes are listed in Supplementary Table 21. A malignant B-cell signature score was calculated in B cells of all nine FL samples using the GSVA package (v.1.38.2) 85 and depicted using the FeaturePlot and VlnPlot functions of Seurat.
Intercellular ligand-receptor interaction analysis. We investigated interactions between NHC subclusters and malignant B cells of nine FL samples (FL 2-FL 10) using the CellPhoneDB package (v.2.1.1) 52 on Python (v.3.6). Gene expression information relevant to each NHC subcluster in integrated FL NHCs was used for NHC data, whereas gene expression information relevant to malignant B-cell clusters in each FL sample was separately used for malignant B-cell data, as gene expression profiles of malignant B cells vary greatly among samples. We then performed pairwise comparisons between NHC subclusters and malignant B-cell clusters. In brief, we derived potential ligand-receptor interactions based on the expression of a receptor gene by one lineage subpopulation and a ligand gene by another. We filtered genes expressed in >20% of cells in any given subpopulation. We then permuted the cluster labels of all input cells 1,000 times and calculated the mean interaction score (the average receptor expression level in a subpopulation multiplied by the average ligand expression level in the interacting subpopulation), which generated a null distribution of the mean interaction score for each ligandreceptor pair in each pairwise comparison across subpopulations. Thereafter, we located observed mean interaction scores that were the same or higher than the actual mean score in the null distribution and calculated the proportion of the observed scores, conferring a P value for the likelihood of specificity of a given ligand-receptor complex to a given cluster pair. To consider interactions between FL NHCs and FL malignant B cells, we selected only interactions with a P value of <0.05 in more than half of FL cases (>4 cases). Furthermore, to assess subcluster-specific lymphomagenesis mechanisms in FL stroma, we extracted interactions that included a molecule in which gene expression was significantly upregulated in at least one FL NHC subcluster compared with that in the corresponding MFLN subcluster. We integrated interaction scores and P values of interactions between pairs consisting of the same NHC subcluster and malignant B-cell clusters from different FL samples, as previously described 84 . In brief, we calculated mean interaction scores for pairs that included the same NHC subcluster and malignant B-cell clusters from different FL samples, then normalized the mean interaction scores per interaction. We also combined P values of interactions for pairs that consisted of the same NHC subcluster and malignant B-cell clusters from different FL samples using Fisher's method. The P values were corrected using the Benjamini-Hochberg method. In Fig. 6a, circles are coloured when gene expression for the indicated stroma-derived factor is upregulated in relevant FL subclusters compared to that in the MFLN counterparts (log fold-change >0 and adjusted P value <0.05).
IF staining. Human LN and lymphoma samples were immediately embedded in OCT compound (Sakura Finetek Japan, 45833) and frozen in hexane cooled with dry ice. Samples were sliced to 3-μm thickness with a cryostat at −20 °C. Sections were dried for 1 h at 20 °C, fixed for 10 min in 4% paraformaldehyde, incubated for 10 min with 0.1% Triton X-100 (Sigma-Aldrich, T9284) for permeabilization, and then treated with 10% goat serum (Sigma-Aldrich, G9023) in PBS or serum-free protein blocking buffer (Dako, X0909) (when using non-goat-derived secondary antibodies) for 30 min. Sections were stained overnight at 4 °C with primary antibodies listed in Supplementary Table 22. After several washes with tris-buffered saline with tween 20 (Sigma-Aldrich, P9416), sections were stained for 1 h with combinations of the following secondary antibodies at 20 °C: AF488-goat-anti-rat IgG (ThermoFisher Scientific), AF594-goat-anti-rabbit IgG (ThermoFisher Scientific), AF594-donkey-anti-goat IgG (ThermoFisher Scientific) and AF647-goat-anti-mouse IgG (ThermoFisher Scientific). A TrueVIEW Autofluorescence Quenching kit (Vector, SP-8500) was used to decrease possible tissue autofluorescence per the manufacturer's instructions. Sections were then mounted in mounting medium with 4,6-diamidino-2-phenylindole (DAPI; Vector, H-1200). Stained samples were imaged using a Leica DMi8 S Platform with the Thunder imaging system (3D Live Cell & 3D Cell Culture & 3D assay). Analysed LNs were verified as malignancy-free by pan-cytokeratin staining. Quantitative analysis of acquired images was performed using ImageJ software (National Institute of Health, v.2.1.0). As LNs and FL carry localized structures, we randomly acquired at least five different regions of interest within each sample and used the median values for statistical analysis.

Flow cytometry analysis of FL haematopoietic cells.
To analyse the expression of CD27 in malignant FL B cells and to perform the binding/adhesion assays described below, we used additionally collected cryopreserved FL samples (FL 11-FL 18). Clinical characteristics of patients in the additional FL cohort are described in Supplementary Table 23. After thawing, cells were filtered through a 70-μm mesh and incubated with PE-anti-CD27 (BioLegend; 1:500), FITC-anti-CD3 (BioLegend; 1:500), APC-anti-CD19 (Miltenyi Biotec; 1:500) and PE-Cy7-anti-CD10 (BioLegend; 1:500) antibodies for 20 min on ice. Cells were then incubated with 7-AAD viability staining solution for 10 min in the dark and analysed using a FACSAria II or III and FlowJo software.
Ex vivo cell adhesion assay. Frozen FL sections were sliced at 6-μm thickness immediately before the assay. For malignant B-cell isolation, we used FL samples in which >90% B cells were confirmed to be malignant by flow cytometry. B cells were isolated from the FL haematopoietic cell suspension using an EasySep Release Human CD19 Positive Selection kit (StemCell Technologies, ST-17754). Cells were then treated with anti-CD27 blocking antibody or isotype mouse IgG1 for 30 min at 4 °C. Thereafter, 2 × 10 6 cells were applied on the sections and incubated with 60 r.p.m. rotation for 5 min, followed by incubation without rotation for 15 min. The incubation with and without rotation was repeated two more times. After incubation, the sections were gently washed with PBS, sealed with a cover glass and imaged using a Keyence BZ-X710 microscope (Keyence). Adherent cells were manually counted using ImageJ.

Prognostic analysis of stroma-derived markers in FL.
To analyse the prognostic potential of gene expression patterns of NHCs in patients with FL, we used a bulk microarray dataset of 180 FL biopsy samples from independent, newly diagnosed cases 66 . To narrow candidates to stroma-specific genes, we initially selected DEGs upregulated in FL BEC and NESC subclusters relative to MFLN counterparts. These were narrowed down to those showing a log fold-change of >0.5 and < 0.1% of cells with an expression level higher than 0 in FL haematopoietic cells (Extended Data Fig. 9a). We did not use genes upregulated in FL LEC subclusters, as the proportion of FL LECs was considerably decreased relative to MFLN LECs and the specificity of these genes to FL stroma was considered unlikely in analyses of bulk tissues. Next, we tested the expression of all candidate genes using the Kaplan-Meier method and two-sided log-rank test. Cut-off expression values of each gene for the Kaplan-Meier survival curves was determined using maximally selected rank statistics 86 . As many putative stroma-specific genes were upregulated in FL, it was possible that P value collection (for example, the Bonferroni method) greatly reduced the number of candidate genes, considering that the sample size in the dataset was not particularly large. Therefore, we extracted genes with reliable prognostic impacts using another approach (Extended Data Fig. 9a,c). We initially divided patients into three groups according to survival outcomes: a favourable group, which comprised patients alive 10 years after diagnosis; an unfavourable group, which comprised patients who died within 5 years of diagnosis; and an intermediate/indefinite group, which comprised the remaining patients (Extended Data Fig. 9c). We then compared the proportion of patients with higher expression of each candidate gene between favourable and unfavourable groups. Genes were considered prognostic when the proportion was significantly higher in the unfavourable group compared with that in the favourable group (Extended Data Fig. 9d). These prognostic genes were further subjected to multivariate analysis (Extended Data Fig. 9a).
To evaluate the prognostic efficiency of the FL TRC signature, we extracted the DEGs that were upregulated in FL TRCs in comparison to MFLN TRCs (Supplementary Table 17). We considered the DEGs with an expression level higher than 0 in <0.1% FL haematopoietic cells, <10% FL BECs and <10% FL LECs and were detectable in the microarray dataset 66 . A total of 11 extracted genes constituting the FL TRC signature are listed in Supplementary Table 19.

nature research | reporting summary
April 2020 Corresponding author(s): Mamiko Sakata-Yanagimoto Last updated by author(s): Jan 20, 2022 Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted Data Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability scRNA-seq data that support the findings of this study have been deposited at the European Genome-Phenome Archive (https://ega-archive.org) database and can be retrieved using the accession number EGAD00001008311. For survival analysis, a DNA microarray dataset from Leich et al66 was downloaded from the Gene Expression Omnibus (GEO) (accession number: GSE16131). For mapping of scRNA-seq data, GRCh38 (https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.39) was used. All other data are available from the corresponding authors on reasonable request. Source data are provided with this paper.

Field-specific reporting
Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences
For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf

Life sciences study design
All studies must disclose on these points even when the disclosure is negative.

Sample size
No statistical method was used to determine sample size a priori. The number of human lymph node and lymphoma samples was highly restricted due to the limited availability of these samples in clinical settings.
Data exclusions Pre-processed single-cell data from each sample were further processed and analysed individually using R package Seurat on RStudio. After removing ribosomal genes, genes expressed in fewer than 3 cells, and cells expressing fewer than 200 genes, we filtered out cells with less than 200 unique feature counts (low quality cells). Cells with unique feature counts greater than three times the median value (possible doublets) and/or cells with more than twice the median number of mitochondrial genes (possible apoptotic or lysed cells) were also removed. After the data integration and clustering analysis, we removed data of NHC subclusters considered possible doublets as characterized by high expressions of marker genes for different NHC components and incongruously high unique feature counts.