Functional heterogeneity of lymphocytic patterns in primary melanoma dissected through single-cell multiplexing

  1. Francesca Maria Bosisio  Is a corresponding author
  2. Asier Antoranz  Is a corresponding author
  3. Yannick van Herck
  4. Maddalena Maria Bolognesi
  5. Lukas Marcelis
  6. Clizia Chinello
  7. Jasper Wouters
  8. Fulvio Magni
  9. Leonidas Alexopoulos
  10. Marguerite Stas
  11. Veerle Boecxstaens
  12. Oliver Bechter
  13. Giorgio Cattoretti
  14. Joost van den Oord
  1. Laboratory of Translational Cell and Tissue Research, KU Leuven, Belgium
  2. Pathology, Department of Medicine & Surgery, Università degli studi di Milano-Bicocca, Italy
  3. ProtATonce Ltd, Greece
  4. National Technical University of Athens, Greece
  5. Laboratory of Experimental Oncology, KU Leuven, Belgium
  6. Department of Surgical Oncology, KU Leuven, Belgium
20 figures, 2 videos, 3 tables and 11 additional files


Single-cell analysis scheme.

(A) TMA construction, Multiplex-stripping immunofluorescence: 60 cores were obtained from 29 patients, assembled in a Tissue Micro Array, and analysed using the MILAN technique; the immunofluorescence images from one round of staining (three markers/round: S100-blue, CD3-red, CD4-green) of three representative cores from our data set are shown. (B) CD8+ cells were analyzed using image analysis for functionality using an activation parameter derived from multiple activation and exhaustion markers evaluated at single-cell resolution; the CD8+ cells are here digitally reconstructed for each of the above standing cores, preserving their spatial distribution in the tissue section and assigning each of them a color according to the activation status. (C) All the cell populations in the cores were phenotypically identified using consensus between three clustering methods and manual annotation from the pathologists; The heatmaps with the levels of expression of the phenotypic markers per cluster for one of the three clustering methods are shown on the left; on the right, all the inflammatory cells are assigned a color based on their phenotype and the tissue is digitally reconstructed for each of the above standing cores, preserving the spatial distribution of each cell. (D) The social network of the cells was analysed using a permutation test for neighborhood analysis in order to make inferences on cell-cell interactions. The results of the neighborhood analysis are generated as a heatmap were the type and the strength of the interaction is expressed with a color code; to simplify the visualization of the interactions, the different cell types are represented in a circle and connected with lines that clarify the type of relationship between them.

Figure 2 with 5 supplements
Definition of activation and implications in overall survival.

A biplot showing the projection of the cells as well as the rotation vectors of the markers over PC2 and PC3 has been created using only CD8-positive cells and the four markers relevant for their functional status: CD69, OX40, LAG3, and TIM3. (A). This was the first step to define a gradient of activation going from the maximum projected value of CD69 (maximum activation) to the maximum projected value of TIM3 (maximum exhaustion) (B). (C) Z-scores of the original markers over PC2 and PC3. (D) Visual representation of the inter- and intra-patient heterogeneity, that shows how most of the patients present a relative homogeneous activation status of the Tcy. Each core is assigned an activation status (‘Active’, ‘Transition’, or ‘Exhausted’). The cores are grouped for each patient, giving an at-a-glance representation of the heterogeneity of the activation status in different areas of the melanoma in the same patient. (E) The survival analysis in our data set shows a higher overall survival for brisk patients (left) and for patients with high levels of activation (right). Most importantly, the functional definition of activation/exhaustion shows improved prognostic value when compared to the brisk morphological parameter (p.value = 0.075 vs p.value = 0.31 log-rank test).

Figure 2—figure supplement 1
Validation of the prognostic implications of our activation score in an independent cohort.

(A) Selection of the data sources: The Cancer Genome Atlas Skin Cutaneous Melanoma (TCGA-SKCM) dataset was deconvoluted using Cibersort with custom cell-specific profiles derived from Tirosh et al [29]. (B) Selection of CD8+ Tcells: First, single cells from Tirosh [29] were mapped (spearman correlation) to the signature profiles included in the LM22 gene signatures from Cibersort. The heatmap represents the contingency matrix between the labels obtained from Tirosh and the ones obtained using the signature profiles of LM22. Cells mapped as T.cells.CD8 were filtered for further analysis. (C) Classification of T.cells.CD8 into ‘Active’ and ‘Exhausted’: our T cell activation/exhaustion scoring system was applied to the CD8+ Tcells classifying them into ‘Active’ and ‘Exhausted’ (see Materials and methods, Functional analysis of TILs). The scatter plot represents the obtained distribution of activation values for each of the cells labeled as T.cells.CD8 (similar to Figure 2.B). (D) Construction of T.cells.CD8 Active and Exhausted gene expression profiles: Custom cell-type-specific profiles were built for T.cells.CD8.Active and T.cells.CD8.Exhausted using the average expression of all the T.cells.CD8 labeled as ‘Active’/‘Exhausted’ (Supplementary file 3). The barplot represents the value for each gene in the profile. (E) Deconvolution of bulk rna-seq data and evaluation of prognosis: TCGA-SKCM dataset was filtered using only stage I and II patients for comparability with our dataset (see Materials and methods, Functional analysis of TILs). Cibersort was applied using the generated profiles obtaining the relative number of T.cells.CD8.Active and T.cells.CD8.Exhausted. Patients were labeled as ‘Active’ if the relative number of T.cells.CD8.Active was higher or equal than the relative number of T.cells.CD8.Exhausted and ‘Exhausted’ otherwise. Survival analysis revealed that the Active group of patients had better prognosis than the exhausted group of patients (log-rank p-value=0.0082) validating the results obtained with our dataset as shown by the Kaplan-Meier curves.

Figure 2—figure supplement 2
qPCR and shotgun proteomics.

Confirmation of the functional subgroups by qPCR (A) and shotgun proteomics plus pathway analysis (B). With qPCR we were able to identify, both in B and NB metastasis, active (green rectangles) and exhausted (red rectangles) cases. With proteomics we confirmed that the ‘IFNg-high’ case had more proteins that the ‘LAG3-high’ and the ‘none’ case involved in inflammatory pathways among which the IFNg-related pathway (in gray).

Figure 2—figure supplement 3
Definition of activation.

The rotation matrix from the PCA shows that PC1 is capturing general expression of the initial markers (all the markers with the same sign) (A). PC1 explained 39,39% of the total variance and was not associated to a patient effect as shown by (B). The 3D animation GIF visible in the on line version of this paper shows how PC1 can be attributed to the general expression of cells some cells express more than others (C): if we visualize a 3D scatter plot with PC1, PC2 and PC3 with every cell colored by activation we can see that the resulting 3D structure has a pyramidal shape with PC1 parallel to its main axis and its apex (main vertex) found at the maximum value of PC1 (Animation 1, Animation 2). As seen in 3D plot and supported by the projections on Figure 2C, points near the vertex (centroid of the projection in Figure 2C), do not show a strong expression for any marker. On the other hand, points close to the base of the pyramid (low values of PC1) show a strong differential expression of the activation and exhaustion markers. Therefore, PC1 was disregarded for the activation/exhaustion analysis and only PC2 and PC3 were used.

Figure 2—figure supplement 4
Definition of Activation in the Core Level (A) and Patient Level (B).

(A) After assessing the activation level of every CD8+ lymphocyte present in the TMA (Materials and methods: functional analysis of TILs and supplementary information 1), a label was assigned to each core (‘Active’, ‘Transition’, and ‘Exhausted’) using one-tailed t-tests comparing the distribution of the activation values in the cells from a particular core versus the background distribution (combination of all cores). p-Values were adjusted using the False Discovery Rate (FDR) method. A cut-off value of 0.001 over the adjusted p-values was used as classification threshold. Histograms show the distribution of cells for a particular core colored in red/dark gray/blue if they were labeled as active/transition/exhausted. (B) In order to obtain patient-specific read-outs, we pooled together all the cells from all the cores for each patient and repeated the same that we did for the cores (to compare the distribution of cells for a specific patient against the background distribution). This classified each patient into one of the three categories (‘Active’, ‘Transition’, and ‘Exhausted’).

Figure 2—figure supplement 5
Biplot showing the projection of the Tcy cells over PCs 2 and 3.

In order to verify the solidity of the method, we compared the results obtained with the DAPI mask with the ones obtained with the CD8+ mask (Figure 2.A). We applied the same approach to define the activation status of the T cytotoxic cells on the Tcy phenotypic cluster, and we observed the same structural behavior of the activation and exhaustion markers (CD69 and TIM3 at opposite ends of the activation spectrum, concordance between LAG3 and OX40).

Activation status and neighborhood analysis in late, early and no regression.

(A) The histogram shows the distribution of the different cores according to the activation levels of the Tcy. The color code identifies the presence and the type of regression areas in the cores. The cases with late regression are all in the left part of the histogram, showing higher levels of activation compared to cores with early regression or without regression. (B) This can be visualized also as a box plot. The neighborhood analysis for late (C) and early (D) regression shows an enrichment in immune-stimulating interactions in the first and more interactions leading to immune impairment in the latter. The thickness of the edge in the network represents the level of interaction between the different cell types. The colour of the line indicates interactions leading to immune suppression (red), to immune stimulation (green), to a probably sub-optimal/impaired immune stimulation (orange), no immune implications (blue).

Figure 4 with 1 supplement
Distribution of the immune cells’ subgroups and significant differences in the morphological and functional TILs categories.

The percentage of cells for each inflammatory subpopulation across all the cores is showed in (A). The 17 phenotypic clusters are on the upper side of the graph, while at the bottom each of them is subdivided in the respective functional subclusters. BC = B cells; cDC1 = classical dendritic cells type 1; cDC2 = classical dendritic cells type 2; Epith = epithelial cells; PC = plasma cells; Lang = Langerhans cells; LV = lymph vessels; Macroph = macrophages; pDC = plasmocytoid dendritic cells; S-M+=S100+MelanA- melanoma cells; S+M-=S100-MelanA+=melanoma cells; S+M+=S100+MelanA+ melanoma cells; Tcy = cytotoxic T cells; Tfh = T follicular helpers; Th = T helpers; Treg = regulatory T cells; suffix: ‘prolif’=proliferating, IFNg = interferon gamma. (B) Significant differences (p.value <0.05) in cell percentages between brisk and non-brisk categories (Wilcoxon rank sum test). (C) Significant differences (p.value <0.05) in cell percentages across the functional groups: Active, Transition, Exhausted (Kruskal-Wallis rank sum test).

Figure 4—figure supplement 1
Phenotype Identification.

A two-tier approach was implemented for the identification of cell subpopulations. (A) Initially, heatmaps generated with KMeans, PhenoGraph, and ClusterX (left) were used to identify the main inflammatory subpopulations, resulting in clusters that were named by manual annotation by the pathologist based on the level of expression of the phenotypic markers; on the right one of the heatmaps is enlarged to serve as an example: for each cluster identified by phenograph the levels of expression of the 15 most expressed markers have been used by the pathologist to identify the inflammatory population most probably present in that cluster. The final phenotypes were defined using a consensus-based approach (a phenotype was assigned if two or more clustering methods agreed in their prediction). (B) Cell populations were further sub-clustered into functional subgroups using PhenoGraph over the set of functional markers and annotation by the pathologist. For each cluster identified with the Materials and method described, here the identified subclusters are shown as an heatmap with the corresponding levels of expression of the selected functional markers.

Neighborhood analysis in the morphological and functional TILs categories.

The result of the neighborhood analysis for active (A), exhausted (B), brisk (C), and non-brisk (D) cases represented as cellular social networks. The thickness of the edge in the network represents the level of interaction between the different cell types. The colour of the line indicates interactions leading to immune suppression (red), to immune stimulation (green), to a probably sub-optimal/impaired immune stimulation (orange), no immune implications (blue). The brisk and active plots display a lot of similairies; nevertheless, the brisk and non brisk categories are not exactly overlapping with the active and exhausted plots, suggesting them to be the result of cases with activation and cases with exhaustion pooled together under the morphological labels. (E) The histogram shows the alternative approach of neighbourhood analysis tailored to explore specifically the interactions regarding melanoma cells. We observed that the main inflammatory cells subtypes in contact with melanoma cells are macrophages and (as expected) epithelial cells, both in brisk and non-brisk cases, followed by Tcy with active and in transition Tcy in brisk cases and proliferating and anergic Tcy in non-brisk cases. Other small differences between brisk and non-brisk cases are more TIM3+ cDC1, cDC2 and TIM3+ macrophages in contact with melanoma cells in brisk cases.

Complete digitalization of a core and practical example of a possible downstream analysis.

Anti-clockwise: (1) images for the markers stained with the MILAN multiplexing technique (for this work 39 markers, but we reached in our laboratory an output of around 100 markers per single section) are acquired and composite images with some selected markers can be prepared. (2) All the markers are used to phenotypically identify all the cell subtypes present in the tissue and the tissue is digitally reconstructed. (3) Studies of the functional status of these cells can be done, for example we could localize exhausted and active cells in the tissue. (4) With neighborhood analysis, it is possible to identify the cells that are in proximity with each other more than chance can suggest, inferring possible interactions between these cell types (in the image, two cells identified as ‘neighbors’ are encircled in red).

Author response image 1
Author response image 2
Author response image 3
Author response image 4
Author response image 5
Author response image 6
Author response image 7
Author response image 8
Author response image 9
Author response image 10
Author response image 11
Author response image 12
Author response image 13
Resampling analysis.

CD8+ TCell populations were sampled 100 times for different sampling sizes (10 to 1000 with a step of 10). For every simulation, we reclassified the patient into “Active”/”Transition”/”Exhausted” using the approach described in the methods. This analysis shows the robustness of the patient classification methodology with most of the patients requiring a very small sample size in order to obtain at least 95% of consistent classifications.

Author response image 14


Animation 1
3D scatter plot with PC1, PC2 and PC3.

Animated in order to see the pyramidal shape of the 3D structure (rotating around the PC2 axis).

Animation 2
3D scatter plot with PC1, PC2 and PC3.

Animated in order to see the pyramidal shape of the 3D structure (rotating around the PC1 axis).


Table 1
Statistical analysis.

Histopathological parameters were correlated with core status (Active/Transition/Exhausted) using pairwise t-tests with pooled standard deviation. Several histopathological parameters were correlated to the average level of activation of each core.

ParameterComparisonStatistical testMultiple testing
p value
Brisk InfiltrationBrisk vs Non Briskt testNo0.8453
Brisk InfiltrationBrisk vs Brisk In Non Briskpairwise t testYes1
Brisk InfiltrationBrisk vs Non Brisk In Non Briskpairwise t testYes1
Brisk InfiltrationBrisk In Non Brisk vs Non Brisk In Non Briskpairwise t testYes1
RegressionRegression vs No Regressiont testNo0.6275
RegressionEarly Regression vs No Regressionpairwise t testYes0.329
RegressionLate Regression vs No Regressionpairwise t testYes0.031*
RegressionLate Regression vs Early Regressionpairwise t testYes0.022*
Count LymphocytesNumber of Lymphocytes vs Level of Activationlinear modelNo0.3714
HistotypeLMM vs NMMpairwise t testYes0.027*
HistotypeLMM vs SSMMpairwise t testYes0.021*
HistotypeNMM vs SSMMpairwise t testYes0.761
Breslow ThicknessBreslow vs Level of Activationlinear modelNo0.9883
UlcerationPositive vs Negativet testNo0.7252
Number of MitosesMore than 6 vs 1 to 6pairwise t testYes1
Number of MitosesMore than 6 vs 0pairwise t testYes1
Number of Mitosesone to 6 vs 0pairwise t testYes1
Key resources table
Reagent type
(species) or
DesignationSource or referenceIdentifiers
AntibodyCD4 rabbit monoclonal EPR6855Abcam1 µg/ml
AntibodyHLA-DR mouse monoclonal IgG2b SPM288AbcamAB_11252171 µg/ml
AntibodyTAP2 mouse IgG1 monoclonal TAP2.17Abcam1 µg/ml
AntibodyCD141 rabbit monoclonal EPR4051AbcamAB_22018051 µg/ml
AntibodyMYC rabbit monoclonal EP121Sigma Aldrich1 µg/ml
AntibodyFOXP3 mouse monoclonal IgG1 236A/E7AbcamAB_4452841 µg/ml
AntibodyMX1 Rabbit polyclonalAbcamAB_106789251 µg/ml
AntibodyLAG3 mouse monoclonal IgG1 11E3AbcamAB_7761021 µg/ml
AntibodyPD-L1 rabbit monoclonal 28–8Abcam/EpitomicsAB_26878781 µg/ml
AntibodyCD1a rabbit monoclonal EP3622Abcam/EpitomicsAB_6269571 µg/ml
AntibodyCD123 mouse monoclonal IgG2b NCL-L-CD123Leica-Microystem/NovocastraAB_105552711 µg/ml
AntibodyPhospho-Stat1 rabbit monoclonal 58D6Cell SignalingAB_5612841 µg/ml
AntibodyCD20 mouse monoclonal IgG2a L26DakoAB_7820241 µg/ml
AntibodyCD1a mouse monoclonal IgG1 O10Dako1 µg/ml
AntibodyCD1c mouse monoclonal IgG1 2D4DakoAB_26230491 µg/ml
AntibodyPRDM1 rat monoclonal 6D3DakoAB_1 µg/ml
AntibodyS100AB rabbit polyclonalDako1 µg/ml
AntibodyCD56 mouse monoclonal IgG1 123C3.D5NeomarkersAB_6271271 µg/ml
AntibodyKi-67 mouse monoclonal IgG2a UMAB107OrigeneAB_26291452 µg/ml
AntibodyLysozyme rabbit polyclonalOrigeneAB_10047661 µg/ml
AntibodyPD-1 mouse monoclonal IgG2a UMAB197OrigeneAB_26291981 µg/ml
AntibodyTIM3 goat polyclonalR and DAB_3552351 µg/ml
AntibodyCXCL13 mouse monoclonal IgG1 53610R and DAB_20860491 µg/ml
AntibodyOX40 mouse monoclonal IgG1 Ber-ACT35Santa CruzAB_6268971 µg/ml
AntibodyIRF4 goat monoclonal M-17Santa CruzAB_21271451 µg/ml
AntibodycMAF rabbit monoclonal M-153Santa CruzAB_6385621 µg/ml
AntibodyBCL6 rabbit monoclonal N3SCBTAB_11580741 µg/ml
AntibodyCD16 mouse monoclonal IgG2a 2H7SCBTAB_5635081 µg/ml
AntibodyCD68 mouse monoclonal IgG3 PGM1Thermo FisherAB_109795581 µg/ml
Antibodyp16 mouse monoclonal IgG2a JC8SCBTAB_7850181 µg/ml
AntibodyMelanA rabbit monoclonal A19-PNovusBioAB_19872851 µg/ml
AntibodyPodoplanin rat monoclonal IgG2a NZ-1.2Sigma AldrichAB_109205771 µg/ml
AntibodyCD69 rabbit polyclonalSigma AldrichAB_26811571 µg/ml
AntibodyCD3 rabbit polyclonalSigma AldrichAB_23356771 µg/ml
AntibodyGBP1 rat monoclonal 4D10Sigma AldrichAB_8289641 µg/ml
AntibodyLangerin rabbit polyclonalSigma AldrichAB_10784531 µg/ml
AntibodyIRF8 rabbit polyclonalSigma AldrichAB_18519041 µg/ml
AntibodyCD8 rabbit monoclonal SP16Thermo FisherAB_6272111 µg/ml
AntibodyCD138 mouse monoclonal IgG1 MI-15Thermo FisherAB_109870191 µg/ml
Table 2
Multiplex antibody panel description.
CD41 µg/mlrabbit MabEPR6855AbcamN/A
HLA-DR1 µg/mlmouse IgG2bSPM288Abcam1125217
TAP21 µg/mlmouse IgG1TAP2.17AbcamN/A
CD1411 µg/mlrabbit MabEPR4051Abcam2201805
MYC1 µg/mlrabbit MabEP121Sigma AldrichN/A
FOXP31 µg/mlmouse IgG1236A/E7Abcam445284
MX11 µg/mlRabbitAbcam10678925
LAG31 µg/mlmouse IgG111E3Abcam776102
PD-L11 µg/mlrabbit Mab28–8Abcam/Epitomics2687878
CD1a1 µg/mlRb mAbEP3622Abcam/Epitomics626957
CD1231 µg/mlmouse IgG2bNCL-L-CD123Leica-Microystem/Novocastra10555271
Phospho-Stat11 µg/mlrabbit Mab58D6Cell Signaling561284
CD201 µg/mlmouse IgG2aL26Dako782024
CD1a1 µg/mlmouse IgG1O10DakoN/A
CD1c1 µg/mlmouse IgG12D4Dako2623049
PRDM11 µg/mlrat6D3Dako628168
S100AB1 µg/mlrabbitDakoN/A
CD561 µg/mlmouse IgG1123C3.D5Neomarkers627127
Ki-672 µg/mlmouse IgG2aUMAB107Origene2629145
Lysozyme1 µg/mlrabbitOrigene1004766
PD-11 µg/mlmouse IgG2aUMAB197Origene2629198
TIM31 µg/mlgoatR and D355235
CXCL131 µg/mlmouse IgG153610R and D2086049
OX401 µg/mlmouse IgG1Ber-ACT35Santa Cruz626897
IRF41 µg/mlgoatM-17Santa Cruz2127145
cMAF1 µg/mlrabbitM-153Santa Cruz638562
BCL61 µg/mlrabbitN3SCBT1158074
CD161 µg/mlmouse IgG2a2H7SCBT563508
CD681 µg/mlmouse IgG3PGM1Thermo Fisher10979558
p161 µg/mlmouse IgG2aJC8SCBT785018
MelanA1 µg/mlrabbit MabA19-PNovusBio1987285
podoplanin1 µg/mlrat IgG2aNZ-1.2Sigma Aldrich10920577
CD691 µg/mlrabbitSigma Aldrich2681157
CD31 µg/mlrabbitSigma Aldrich2335677
GBP11 µg/mlrat4D10Sigma Aldrich828964
Langerin1 µg/mlrabbitSigma Aldrich1078453
IRF81 µg/mlrabbitSigma Aldrich1851904
CD81 µg/mlrabbit MabSP16Thermo Fisher627211
CD1381 µg/mlmouse IgG1MI-15Thermo Fisher10987019

Additional files

Source code 1

Code for all the analyses included in the paper.
Source data 1

Image metadata.
Source data 2

Annotated single cell data and patient survival.
Source data 3

Correlation histopathological data and lymphocyte features.
Supplementary file 1

Contingency table between the cell labels given by Tirosh et al. (2016) and the predicted cell type from the LM22 signatures from Cibersort.

The highest spearman correlation was obtained for the predictions.
Supplementary file 2

Transcription factors selected to represent the proteins used to calculate the activation score.
Supplementary file 3

Gene expression signatures for the T.cells.CD8.Active and T.cells.CD8.Exhausted cell types.
Supplementary file 4

Relative percentages of T.cells.CD8.Active and T.cells.CD8.Exhausted cell types in the TCGA-SKCM dataset (stages I and II).
Supplementary file 5

Justification of the statistical test chosen for each analysis.
Supplementary file 6

Patient metadata.
Transparent reporting form

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Francesca Maria Bosisio
  2. Asier Antoranz
  3. Yannick van Herck
  4. Maddalena Maria Bolognesi
  5. Lukas Marcelis
  6. Clizia Chinello
  7. Jasper Wouters
  8. Fulvio Magni
  9. Leonidas Alexopoulos
  10. Marguerite Stas
  11. Veerle Boecxstaens
  12. Oliver Bechter
  13. Giorgio Cattoretti
  14. Joost van den Oord
Functional heterogeneity of lymphocytic patterns in primary melanoma dissected through single-cell multiplexing
eLife 9:e53008.