Measuring Intratumoral Heterogeneity of Immune Repertoires

There is considerable clinical and fundamental value in measuring the clonal heterogeneity of T and B cell expansions in tumors and tumor-associated lymphoid structures—along with the associated heterogeneity of the tumor neoantigen landscape—but such analyses remain challenging to perform. Here, we propose a straightforward approach to analyze the heterogeneity of immune repertoires between different tissue sections in a quantitative and controlled way, based on a beta-binomial noise model trained on control replicates obtained at the level of single-cell suspensions. This approach allows to identify local clonal expansions with high accuracy. We reveal in situ proliferation of clonal T cells in a mouse model of melanoma, and analyze heterogeneity of immunoglobulin repertoires between sections of a metastatically-infiltrated lymph node in human melanoma and primary human colon tumor. On the latter example, we demonstrate the importance of training the noise model on datasets with depth and content that is comparable to the samples being studied. Altogether, we describe here the crucial basic instrumentarium needed to facilitate proper experimental setup planning in the rapidly evolving field of intratumoral immune repertoires, from the wet lab to bioinformatics analysis.


INTRODUCTION
Analysis of T and B cell repertoires has become a valuable and powerful tool for characterizing the immune response, complementing other high-content approaches such as transcriptome analysis and mass cytometry (1)(2)(3). Both T and B cell repertoires have been shown to be predictive of survival for cancer patients (4)(5)(6)(7)(8) and for assessing response to checkpoint immunotherapy (9)(10)(11)(12)(13).
Immune repertoires allow for tracking lymphocyte lineages, and make it possible to trace the evolution and heterogeneity of anti-tumor immunity. For example, immune repertoire analysis has been used to reveal that most tumor-infiltrating regulatory T cells in human cancers originate in the thymus rather than from local conversion of conventional T cells (14)(15)(16). Tumeh et al. have shown that pre-existing tumor-infiltrating CD8 + T cell clones expand after PD-1 checkpoint blockade and this underlies the positive response of advanced melanoma patients to therapy (10). In contrast, more recent study of basal and squamous cell carcinoma treated with anti-PD-1 revealed clonal expansion of novel clonotypes that had not previously been observed in the tumor, but not of pre-existing tumor-infiltrating lymphocytes (TILs) (17). Based on limited available information on T cell repertoires, Zhao et al. have suggested that there is increased clonality (i.e., decreased diversity) of TILs after anti-PD-1 therapy in treatment-responsive cases of glioblastoma and decreased clonality (i.e., increased diversity) in non-responders (18). On the other hand, Schalper et al., made the determination based on (also limited) TCR repertoire information obtained using multiplex PCR from formalin-fixed paraffin-embedded (FFPE) block-extracted RNA that a decrease in T cell clonality in glioblastoma patients receiving the same treatment was associated with longer survival (19).
In one of the first deep sequencing-based efforts to estimate intratumoral T cell receptor (TCR) repertoire heterogeneity, Emerson and colleagues reported relatively high TCR repertoire similarity throughout ovarian carcinoma tumors in terms of clonal overlap (20). At the same time, in melanoma, it has been shown that TILs harvested from different tumor fragments possess different reactivity against melanoma-associated antigens, and that corresponding epitopes were typically found in the same tumor fragments as their cognate TCRs (21,22), although these clones were present at low frequencies (21). Subsequent studies have revealed intratumoral heterogeneity of T and B cells in lung adenocarcinoma (6,23,24), esophageal squamous cell carcinoma (25,26), colorectal (27), ovarian (28), breast (29), and pancreatic cancers (30).
In several recent studies, immune repertoires were analyzed in conjunction with genomic heterogeneity of tumor cells. These data revealed that immune surveillance evolves with the tumor, with certain T cell clones tracking tumor neoantigens both spatially and temporally and imposing selection pressure on the tumor (6,23,24,27,28). For B cells, the affinity maturation of immunoglobulins against tumor antigens was deduced from repertoire analysis (5,7,31), although no association with tumor clonality and evolution has been revealed to date.
One of the main complications in comparing immune repertoires between different timepoints and distinct tumor sections is that there is always some dispersion in clonal frequencies caused by sampling limitations. These can arise at the level of tissue sampling, T or B cell counts, the amount and quality of extracted genomic DNA (gDNA) or RNA, efficiency of cDNA synthesis and template molecules entrance into PCR amplification. The stochastic nature of PCR amplification further adds artificial dispersion. RNA-based TCR and immunoglobulin repertoire profiling is particularly informative in terms of assessing the functional activity of infiltrating T and B cells and local antibody production, especially given the ultra-high immunoglobulin expression levels seen in plasma cells. However, in such experiments, the cell-to-cell dispersion in expression levels can further increase the artificial repertoire heterogeneity originating from one or more randomly-sampled or undersampled plasma cells or active effector T cells. These sources of natural and technical noise have to be taken into account and distinguished from actual repertoire differences.
Some of the aforementioned works employed replicates produced at the level of split gDNA samples. This potentially allowed the authors to control for the dispersion in TCR or BCR repertoire content arising from technical errors (20,25,26,30), with the exception of sample-to-sample heterogeneities associated with the DNA extraction procedure and thus, importantly, differences in sampling depth. However, this information has not been implemented to build an appropriate noise model.
Here, we propose a straightforward approach for measuring the extent of TCR and immunoglobulin heterogeneity in tumor samples relative to internal controls in terms of baseline repertoire dispersion. This is estimated using tissue sample replicates that have been split at the level of homogenized cells. These replicates are then used to train a beta-binomial noise model, as suggested by Rytlewski et al., which is further used to exclude false positives amongst differentially represented clones (32).
We use this approach to demonstrate uneven clonal distribution of CD8 + T cells in poorly-infiltrated B16F0derived tumors in a mouse model of melanoma (33). This is associated with clustered distribution of these cells within the tumor, based on multicolor fluorescent immunohistochemistry analysis, indicating local proliferation of CD8 + T cells in situ at the tumor site. Working with human tissue samples, we also demonstrate heterogeneous distribution of plasma cell clones in a lymph node heavily infiltrated by metastatic melanoma and in a primary colorectal tumor. We also show a scenario in which training with high quality, deeplyanalyzed biological replicates may lead to identification of false-positive clonal expansions when analyzing more noisy samples of interest. This highlights the importance of replicas for correct repertoire comparison, and of careful use of this analytical tool.

Lymphocyte Infiltration Pattern of B16F0 Melanoma
The spatial clonal heterogeneity of TILs has not been thoroughly studied in mouse tumor models, and it is an intriguing question whether such heterogeneity exists and how it can affect repertoire-based analysis. Uncovering such heterogeneity could also shed light on sources of TILs for corresponding models. In order to reveal possible sources of TIL clonal heterogeneity within tumors, we first studied their patterns of distribution in mouse melanoma. Using multicolor IHC, we analyzed the distribution of CD4 + /CD8 + T cells and B cells in whole tumor tissue slices from a B16F0 melanoma model. We found a common distribution pattern for all lymphocyte subsets, with prominent accumulation in the fibrous tumor capsule and in several large clusters within tumor nodes (Figure 1). The tumor capsule is characterized by a high density of immature, hyper-permeable blood vessels that facilitate lymphocyte infiltration (34), while surrounding loose connective tissue offers a perfect substrate for further lymphocyte migration (35). This may result in relatively non-specific lymphocyte accumulation in the surrounding tumor envelope. Prior work has also shown that T cells in tumor nodes are more clonal and associated with lower clonal diversity compared to stromal T cells in ovarian tumors (28).
Lymphocyte clusters within the tumor were also related to certain morphological structures, as revealed by comparison of fluorescently-labeled and histological slices. One common feature of these structures was the presence of a blood vessel within the pocket that almost exclusively contained leukocytes ( Figure 1C). It should be noted that only about 25% of blood vessels within the tumor were so prominently surrounded by leukocytes. These are likely to be high endothelial venule pockets that have analogous histological appearance, and give rise to tertiary lymphoid structures (36)(37)(38). These intratumoral clusters of CD4 + and CD8 + T cells may originate from locally enhanced infiltration and/or local proliferation of clonal T cell populations. The latter would be expected to lead to a highly heterogeneous distribution of T cell clones across the tumor.

Pipeline for Measuring Heterogeneity and Local T Cell Expansion
To clarify the origin of observed clusters, we designed a pipeline that allows to measure the contribution of local clonal expansions to repertoire heterogeneity. This approach fully accounts for natural dispersion in clonal frequencies between samples that originates from sampling limitations and is unrelated to real clonal heterogeneity (Figure 2). We processed tumor masses using two alternative methods, both of which produce multiple single-cell suspensions. In the first setup (experimental), the tumor is dissected into four fragments of comparable size, and each fragment is processed independently. In the second setup (control), the whole tumor sample is first homogenized and filtered, after which the four replicates are split from the resulting PBSwashed single cell suspension. Further TCR repertoire profiling of these control samples allows to measure the natural sample-tosample variation resulting from stochastic factors associated with cell sampling and sorting, cell-to-cell variation in TCR mRNA expression, and mRNA and cDNA sampling in the course of library preparation.
In these experiments, we were particularly interested in the nature of intratumorally-observed clusters of T cells. Therefore, we carefully cleaned the excised tumors from the collagenous envelope, including the fibrous tumor capsule shown in Figure 1, in order to focus our analysis on lymphocyte heterogeneity within the tumor parenchyma. It should be noted that our ability to estimate clonal frequencies-and thus the extent of correlation of these frequencies between replicates and distinct tumor sections-is intrinsically limited by sampling bottlenecks and the resulting depth of profiling in terms of cell counts and template gDNA/cDNA counts. Table 1 summarizes the number of sorted cells and unique UMI-labeled TCRβ cDNA molecules analyzed for each sample.

Measuring Intratumoral Clonal Heterogeneity of T Cells
Seven melanoma-bearing mice were randomly subdivided into two groups with average tumor volume of 0.23 ± 0.7 (Mean ± SD) mm 3 and 0.32 ± 0.11 mm 3 for control and experimental group, respectively, and processed by either experimental procedure. After removal of outer tumor capsule density of T cell infiltration did not correlate with tumor volume and constituted 280 ± 150 cells per µl of tumor tissue.
We compared the correlation of clonal TCRβ frequencies for the replicates obtained using both setups. As expected, clonal frequencies were highly correlated between biological replicates in control tumor sections ( Figure   most abundant clones were of comparable frequency in all replicates. At the same time, this correlation was not ideal, and the extent of this disparity reflected the natural dispersion between hypothetically "identical" tumor replicates due to sampling effects and stochasticity during library preparation. In contrast, correlation of clonal frequencies between experimental tumor sections was significantly lower (Mann Whitney U-test, p < 0.0001; Figures 3B,C). This poorer correlation compared to the control method reflects the true heterogeneity between the analyzed samples, accounting for all of the technical limitations and bottlenecks. In some of the experimental tumor sections, particular CD8 + clones showed drastic expansion ( Figure 3B). For example, the clonal TCRβ CDR3 variant CGARDWEDAEQFF occupied ∼19% of the repertoire in tumor section #3 but <3% in the other sections. Similarly, CDR3 variant CASGDALGYEQYF occupied ∼13% of the repertoire in section #3 vs. just 0.1-1.8% in the other sections. To statistically identify clonotypes that are differentially expanded across two tumor sections, we used an approach suggested by Rytlewski et al., in which a pair of control repertoires with medium numbers of clones was used to train a beta-binomial noise model. Using this model, we identified significantly contracted and expanded clones between each pair of repertoires ( Figure 3C; Mann-Whitney U-test, p = 0.0001). Notably, such clones were only found in experimental samples, and never between control replicate samples. These results were stably reproducible in three control and four experimental tumor preparations (Figures 3C,D). Based on these findings, we concluded that the observed intratumoral clusters of CD8 + T cells (Figure 1) are probably formed by progeny of T cells that infiltrated the tumor via high endothelial venules, but then proliferated in situ to form relatively large local clonal expansions.
It must also be noted that such local clonal expansions were consistently identified almost exclusively among CD8 + but not among CD4 + T cells. In general, CD4 + T cells rarely form expansions as large as those formed by CD8 + T cells, reflected by lower clonality score (see Table 1). Thus, we believe that the sensitivity and accuracy of repertoire heterogeneity analysis is insufficient to reliably assess relatively minor clonal expansions that would be expected to intrinsically occur amongst CD4 + T cells. Analysis of larger lymphoid structures [e.g., tertiary lymphoid structures in human cancers (38)] should reveal statistically significant local CD4 + expansions.

Measuring Intratumoral Clonal Heterogeneity of B Cells
Heterogeneity of immunoglobulin repertoires across tumor tissues can be investigated in a similar fashion, with the caveat that immunoglobulin expression levels differ dramatically between naive, memory, and plasma B cells by an average ratio of ∼2:5:500 (40), with high dispersion between B cell clones and individual cells. The RNA-based analysis of immunoglobulin repertoires mainly provides functional information on the relative abundance of locally produced clonal antibody variants and the most activated effector B cell receptor (BCR) distributions-mixed together, if functional B cell subsets were not sorted first. Thus, questions related to clonal composition of infiltrating B cells-irrespective of their functional activity-should be preferably studied with DNAbased immunoglobulin profiling.
Here, we studied -at the RNA level -immunoglobulin heterogeneity in a sample from the metastatically-infiltrated lymph node of a patient with cutaneous melanoma. We used a modified experimental scheme that compares two tissue sections, with the advantage of having two internal controls and four independent experimental comparisons (Figure 4). According to flow cytometry analysis, CD45 + leukocytes constituted 17.7 and 4.8% of all cells in analyzed fragments, of which 28 and 25% were CD19 + B cells, respectively. CD20-CD38 + plasma cells constituted 11.5 and 21.6% of all B-cells. Fragments were comparable in size, and the number of cells we isolated was also comparable (7,500 and 4,700 plasma cells in total).
As expected, clonal frequencies correlated well between replicates from the same section ( Figure 5A). In contrast, samples from different tumor fragments showed much lower correlation, as we observed for TCR repertoires from murine tumors. For example, the IGH CDR3 clonal variant CARSGGYFDWGFFDYW occupied 9.9% of the repertoire in tumor fragment X, but represents <0.1% in fragment Y. Likewise, clonal variant CARVGTGTKSFDYW occupied 9% of the repertoire in fragment Y and only 3.5% in fragment X (Figure 5B). Clonal expansions were identified only when comparing the samples obtained from X and Y fragments (Figure 5C), and lower correlation of clonal frequencies between X and Y fragments compared to replicate samples allowed to estimate the level of immunoglobulin heterogeneity ( Figure 5D). Along with the composition, clonality, and hypermutation of intratumorally-produced antibodies, the proportion of distinct antibody isotypes may also be a crucial parameter with prominent prognostic value, as has been shown for human melanoma and subtypes of lung adenocarcinoma and bladder cancer (5,8). Like clonal frequencies, isotype proportions are heterogeneous across tumor tissue and subject to sampling noise, but also generally correlate well between replicates produced at the level of single-cell suspensions ( Figure 5E). Thus, the findings above are equally applicable to the study of the heterogeneity of the isotypic composition of antibodies produced in tumor tissues and the correlation of this heterogeneity with unevenness of the immune landscape.

Noise Models Must Be Trained on Datasets of Comparable Depth
We next repeated the whole pipeline to assess immunoglobulin heterogeneity in human colon cancer sections, using the same experimental design shown in Figure 4. According to flow cytometry estimates CD45 + leukocytes constituted 35-45% of all cells in analyzed section and consisted predominantly from lymphocytes (70-80%), that in turn contained 20-40% of CD19 + B cells. Remarkably, both sections contained unusually high proportions of CD19 + CD20 − CD38 + plasma B cells among all CD19 + B cells, on the order of ∼50%, suggesting the presence of tertiary lymphoid structures (38). It is also important to note here that in spite of similar volume of the two sections, the total number of isolated cells was significantly different (∼1 million vs. 10 million). As a result the two sections contained very different numbers of CD19 + B cells and plasma cells, with about 50,000 in fragment X and about 500,000 in fragment Y ( Table 2). Both replicates from fragment Y contained high numbers of plasma cells, providing excellent statistics for the accurate identification of clonal frequencies (Figure 6A). We trained the noise model on these replicates, and used it to estimate heterogeneity between the X and Y fragments. This analysis identified considerable heterogeneity and lots of clonal expansion, as expected ( Figure 6B). However, when we applied the trained model to the replicates obtained from fragment X, we were disappointed to observe a number of false-positive clonal expansions in both X1 and X2 (Figures 6C,D). Thus, we concluded that since the X1 and X2 replicates were less rich in cells and correlated poorly with the beta-binomial noise model trained on a pair of extra-deep control repertoires (Y1 and Y2), this model must be erroneously identifying statistically significant clonal expansions.
This experiment clearly shows how easy it is to make a mistake in interpreting data pertaining to local immune repertoire clonality, even in an apparently well-controlled experimental setup. Training on replicates of appropriate and comparable depth, both in terms of cell counts and TCR/immunoglobulin cDNA molecule counts, is therefore critical. For identification of immunoglobulin clonal expansions at the mRNA level, which mainly reflect clonality of locally-produced antibodies (5), we recommend training noise models on replicates containing comparable numbers of plasma cells among all samples of interest.
Nevertheless, we observed a prominent difference between the X and Y subsections, and note that the extent of miscorrelation differed significantly between experimental (X1/Y1, X2/Y1, X2/Y1, X2/Y2) but not control (X1/X2, Y1/Y2) samples ( Figure 6E). Thus, the proposed approach is relatively stable against additional noise resulting from limited sampling depth when used to estimate the heterogeneity of immune repertoires between tumor sections based on general miscorrelation of clonal frequencies.

DISCUSSION
It is difficult to measure immune repertoire heterogeneity in terms of relative overlap between repertoires sampled from different parts of the same tumor. Take the example of two tissue sections that initially contain identical repertoires of 100 different T cells where each clone is represented by a single cell. If, in our experimental setup, we were to sample a repertoire of 10 T cells from each section, the overlap would typically comprise 1 TCR variant-i.e., 10% of each sample. If we sampled 50 T cells from each sample, the overlap would grow to 25 variants, or 50% of the repertoire. If we sampled all T cells, the overlap would reach 100%. This simplistic example clearly shows that measurements of repertoire overlaps between tissue sections, T or B cell subpopulations, or time points strongly depend on the depth of repertoire analysis in terms of the number of sampled cells and TCR/immunoglobulin molecules (41). One possible way to analyze relative repertoire overlap across sections is normalization of profiling depth, by downsampling to the same number of analyzed cells or UMI-labeled template TCR or immunoglobulin cDNA/gDNA (39).
Furthermore, a clone present in multiple samples should not be simply defined as being ubiquitously present, since its frequency may differ in terms of order of magnitude. On the other hand, no clone can be defined as non-ubiquitously present, since any given clone could be easily lost in analysis due to sampling limitations.
Prominent cell-to-cell differences in TCR and (especially) immunoglobulin expression levels, inaccuracies in RNA extraction and library preparation, and the stochastic nature of PCR amplification all contribute to degradation of repertoire data quality, and this noise contribution is increased at limited sampling volumes. These considerations dictate that appropriate noise models should be trained on replicates obtained at the level of single-cell suspensions, and that such training should be performed on replicate samples in which T or B cell numbers and expression levels are comparable with those of the samples of interest. The approach that we propose here addresses both issues, enabling (1) estimation of repertoire heterogeneity across and between tumor tissues, and (2) identification and quantification of locally expanded TCR or immunoglobulin clonotypes. The approach controls for the sampling limitations in a given experimental setup, and allows one to quantitatively judge the extent of repertoire heterogeneity in terms of miscorrelation in clonal frequencies between two samples, as compared to the correlation between control replicates obtained at the level of split single-cell suspensions from the same sample.
This approach can also be used to unequivocally identify the number of locally expanded clonotypes and measure the extent of their expansion. Although one cannot conclude whether a particular T or B cell clone is absent in a tumor section under the proposed conditions, it is possible to determine whether a particular clone is locally expanded to a statistically significant extent. The number and size of such local expansions can be measured and compared between different tumor sections, time points, patients, or tumor subtypes.
Linking tumor (42) and immune heterogeneity is critical since local lymphocyte expansions may correlate with the presence or absence of corresponding immunogenic epitopes, that, in turn, may determine the efficiency or inefficiency of immune response against the whole heterogeneous tumor. Therefore, reliable detection of locally expanded T or B cell clones may have important practical applications.
On the one hand, in contrast to the analysis of averaged bulk repertoire of tumor-infiltrating lymphocytes, identification of locally dominant T or B cell clonal expansions may reveal efficient immune response to particular tumor-associated antigens or neoantigens present in a tumor section (21)(22)(23)28). Therefore, reliable capturing of such local expansions could help to reveal tumor-specific T and B cell clones, facilitating development of adoptive T cell and CAR-T therapeutic approaches. Interestingly, detection of tumor-specific TCRs can be further improved if convergent clones targeting the same epitope are found (43).
Second, T and B cell clonal heterogeneity may reflect the overall heterogeneity of immunogenic targets across the tumor. The T cell pressure can sculpt the antigenicity of tumors escaping from immune control (44). At the same time, tumor heterogeneity may be associated with poor prognosis (45,46), either due to its higher evolutionary flexibility or antigen "dilution". The rational way to cope with this intrinsic heterogeneity of permanently evolving tumor is to simultaneously target multiple antigens which distribution across the tumor tissues is non-uniform. Identification and engagement of multiple locally expanded T cell clones from distinct tumor sections could potentially assist in such work.
Along with T cell clonality, B cell expansion was also shown to correlate with response to checkpoint immunotherapy (47,48). Screening of antibodies reconstituted from circulating plasmablasts of responding patients revealed that many of them bind to non-autologous tumor tissues (11). Hence using B cell clones that are locally expanded before or after therapy could substantially narrow the panel of antibodies with potential reactivity against tumor antigens.
The hopes of today's oncology researchers are, to a great extent, connected with progress in the study of immune repertoires, and the development of methods for the reliable and rapid identification of predictive immune signatures and therapeutically relevant T and B lymphocyte clones. The ability to reliably judge the degree of heterogeneity of immune repertoires and capture local clonal expansions will be an important component of these efforts, and we thus hope that our work will become a useful advance along this yellow brick road.

Murine Tumor Model
Experiments were carried out on C57BL/6 female mice aged 3-5 months. Tumors were generated by subcutaneous (s.c.) injection of 10 5 B16F0 cancer cells in 300 µL PBS into the left flank. B16F0 melanoma cells were grown in DMEM medium supplemented with 10% fetal bovine serum (FBS), 0.06% Lglutamine, 50 units/mL penicillin and 50 µg/mL streptomycin. Cells were incubated at 37 • C and 5% CO 2 , and passaged 2-3 times per week. Right before injection, cells were detached by trypsin, counted, and resuspended at a final concentration of 10 6 cell in 3 mL PBS. After 3 weeks, ∼60% of tumors reached a size of ∼1 cm. Mice with linear tumor size ranging from 0.5 to 1.2 cm were sacrificed with isoflurane (Esteve, Italy) in a single day and tumors were surgically removed and prepared for further analysis. All experiments on animals were carried out in accordance with the National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals (NIH Publications No. 8023, revised 1978). The experimental protocol was approved by the Ethical Committee of the Privolzhsky Research Medical University Academy, Russia (EC #6, granted April 17, 2019).

Mouse Melanoma Resection and Lymphocyte Isolation
Freshly-excised tumors were thoroughly cleaned from the outer tumor capsule and either cut into pieces or processed as a whole. For lymphocyte isolation, excised tumor nodules or tumor fragments were homogenized with a gentleMACS dissociator (Miltenyi Biotec, Germany) and incubated in 1-2 mL RPMI medium supplemented with 417 µg/mL Liberase TL (Roche, Germany) and 10 µg/mL DNase I (Roche, Germany) for 30 min at 37 • C in a shaker. After dissociation, cell suspensions were passed through a 70 µm cell strainer and washed twice with 5 mL of incubation buffer (PBS pH 7.2, containing 0.5% bovine serum albumin and 2 mM EDTA). Cell pellets were resuspended in 100 µL of incubation

Immunohistochemical Staining and Analysis
Immunohistochemical (IHC) analysis of mice tumors was done using Zn-fixed paraffin-embedded tissue slices prepared as described previously (49). Briefly, a 2-3-mm-thick piece was cut from the side of the tumor with the cutting plane perpendicular to the skin, and then transferred into formalin-free Zn fixative For IHC staining, slides were dried and marked with a hydrophobic barrier pen (Diagnostic BioSystems), then washed in TBST buffer (Cell Marque, USA) for 5 min and blocked with 5 mg/mL casein (PerkinElmer) in TBST for 30 min at RT. Predissolved primary and secondary antibodies or streptavidin were added directly into blocking solution in a dropwise manner. The following primary antibodies, dilutions, and incubation conditions were used: 1:500 anti-CD4 clone EPR19514 (Abcam), overnight at 4 • C; 1:250 biotinylated anti-CD8 clone 53.6-7 (BioLegend), 2 h at RT; and 1:250 biotinylated anti-CD45R/B220 clone RA3-6B2 (BioLegend), 2 h at RT. After staining, slides were washed twice in TBST for 5 min, blocked with 5 mg/mL casein in TBST for 10 min at RT, and then incubated for 1 h at RT with 1:1,000 HRP-labeled secondary antibodies (PerkinElmer) or 1:200 HRP-labeled streptavidin (PerkinElmer) in blocking solution. After secondary antibodies the slides were then washed twice in TBST for 5 min and incubated with a 1:75 dilution of Opal TSA reagent (PerkinElmer) in amplification diluent (PerkinElmer) for 10 min at RT. Primary antibodies and Opal reagents were used in the following order and combinations: anti-CD4 with Opal520, anti-CD8 with Opal570, and anti-CD45R/B220 with Opal650. Opal-treated slides were washed once with TBST for 5 min and processed for antibody stripping and re-staining. For antibody stripping, slices were submerged for 10 min in 0.1 M glycine solution (pH 10.0) and then washed once in TBST before re-staining. After antibody staining was complete, 3 µg/mL DAPI was applied in PBS buffer for 10 min, followed by a 5 min wash in ddH 2 O. Slices were immersed in glycerol and covered with a coverslip that was fixed with a nail polish. Fluorescent images of whole slices were acquired with a 10X objective (NA0.45) on an Eclipse Ti microscope (Nikon, Japan) with an Andor Neo high dynamic range sCMOS camera (ANDOR, UK). The following excitation and emission filters were employed: 377/50 and 447/60 for DAPI, 480/40 and 535/50 for Opal520, 525/50 and 600/50 for Opal570, and 620/60 and 705/72 for Opal650. Whole slices were scanned in multi-point acquisition mode with NIS Elements software (Nikon, Japan) with 10-20% image overlap. To obtain multicolor images of the whole slice, single-color images were stitched with the MIST plugin for ImageJ (50) and overlaid.
H&E staining was done with Mayer's hematoxylin and eosin (Biovitrum, Russia). Whole tissue slices were scanned with a DM2500 microscope (Leica) equipped with motorized stage (Märzhäuser Wetzlar, Germany) using a 20X objective (NA0.4) and LAS software. Images were stitched with the MIST plugin. For analysis of BCR heterogeneity in human metastatic melanoma, an excised lymph node that was fully invaded by tumor metastasis was first cleaned of surrounding connective and fatty tissue. The tumor was cut into two pieces, which were then processed separately. tumor pieces were cut into smaller pieces and incubated in 1-2 mL RPMI medium supplemented with 417 µg/mL Liberase TL (Roche, Germany) and 10 µg/mL DNase I (Roche, Germany) for 30 min at 37 • C in 5% CO 2 . Samples were then homogenized with a gentleMACS dissociator (Miltenyi Biotec, Germany). After dissociation, cell suspensions were passed through a 70 µm cell strainer and washed twice with 50 mL of PBS (pH 7.2) containing 0.5% bovine serum albumin and 2 mM EDTA. The pellets were resuspended in PBS at a concentration of 5 × 10 5 cells/100 µL. Two replicates of 50 µl cell suspension was lysed in 350 µl of RLT cell lysis. Colorectal tumor sections were processed similarly with minor modifications (e.g., homogenization with gentleMACS was performed before incubation in Liberase TL and DNase I). Two 500-1,000 µm 3 colorectal tumor sections were taken from distant regions of the primary tumor, which were ∼3 cm apart.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found on Figshare (https://figshare.com/s/98d3d72668acabe91a64).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the local ethical committee. The patients/participants provided their written informed consent to participate in this study. The animal study was reviewed and approved by the Ethical Committee of the Privolzhsky Research Medical University Academy, Russia (EC #6, granted April 17, 2019).