From specimen excisions taken for confirmation of the diagnosis or for follow-up of stage ≥IIA MF at the dermatology departments of the University Hospital Essen and the Ruhr-University Bochum, part of the tissue from 12 patients could also be used fresh for further scientific investigations. The ethics committee of the University Duisburg-Essen approved the project (18–8230-BO), which was conducted in accordance with the Declaration of Helsinki. Written informed consent from each patient was obtained prior to any measures.

Native tissue of cutaneous MF lesions was immediately (i.e., within 2 hours after excision) transferred to single-cell suspension. Briefly, tissue cubes of 2 mm edge length were minced in gentleMACS (Miltenyi Biotec, Bergisch Gladbach, Germany) using program “h_tumor_01”, followed by program “h_tumor_02” twice with 30 min intervals in an enzyme mix consisting of 200 µl Enzyme H, 100 µl Enzyme R and 25 µl Enzyme A (catalogue #130–095-929 Miltenyi Biotec). After three washes with PBS/0.05% bovine serum albumin (BSA), cells were passed through a 100 µm cell strainer. The number and viability of T cells in the single-cell suspension was assessed by flow cytometry using antibodies against CD3 (catalogue #1P-514-T025, EXBIO Praha, Vestec, Czech Republic) and CD45 (catalogue #1A-222-T100, EXBIO Praha) together with 7-aminoactinomycin (7AAD, catalogue # SML1633–1ML, Merck, Darmstadt, Germany) in a CytoFLEX (Beckman Coulter, Brea, CA, USA). More than 80% living cells containing at least 30% T cells were required for scRNA and scTCRseq using the Single Cell 5’ Library & Gel Bead Kit v1.1 (catalogue #1000165, 10xGenomics, Pleasanton, CA, USA) for 12,000 cells at a concentration of 600 cells/µl. After Gel Bead-In EMulsions (GEM) generation, the mRNA was reversely transcribed and cDNA amplified for library construction. The gene expression libraries were sequenced at the DKFZ Genomics and Proteomics core facility on a NovaSeq 6000, paired-end with 26 cycles of the forward and 74 cycles on the reverse read, the TCR libraries on a NextSeq 550 platform with Paired-End 150bp Mid-Output. Samples were sequenced at an average of 64,251 reads per cell. The number of recovered cells varied between 1,092 and 9,846 cells with an average of 5,246 cells (43,7% recovery) and 1,480 identified genes per cell.

Sequencing data was processed with the Cell Ranger pipeline (v3.0, 10xGenomics) with default settings respective for gene expression using the hg19 reference genome; TCR libraries were handled accordingly. Quality metrices for filtering living cells were as follows: percentage of mitochondrial gene expression <12%, unique molecular identifier (UMI) count >500, gene count >200 and number of expressed housekeeping genes >40. Number of expressed housekeeping genes was determined by housekeepers identified from Tirosh et al. with at least one UMI count (10). For the subsequent analysis, the Seurat R package (v3.2.2) was used (11). Additional analyses were performed in R (v3.40) using the indicated packages. Samples were individually processed. After log-normalization of gene-expression counts cell cycle effects were regressed using S- and G2M scores depending on markers from Tirosh et al. as well as mitochondrial gene expression (10, 12). The 2000 most variable genes were selected and principal components (PC) were calculated. Between 15 and 25 PCs were used for downstream analysis determined by the inflection point in an elbow plot ranking the PCs by their variance. Cell clusters were determined using shared nearest neighbor modularity optimization using the Louvain algorithm (at resolution of 0.5) and visualized in reduced dimensionality using uniform manifold approximation and projection (UMAP). Significantly differentially expressed genes were assessed using Wilcoxon rank sum test (ln fold change (FC) > |0.25|, Bonferroni-adjusted P value < 0.05). Finally, clonotype information was added to the meta data of the respective Seurat object based on the “filtered_contig_annotations” file. Of note, the identified malignant T-cell clone in some patients was split into different clusters based on the detection of either only α, only β or α and β chains, which were brought together again for subsequent analyses. After processing each sample individually, they were merged with the “merge” command of Seurat and subsequently log-normalized again.

Large-scale chromosomal copy number variations were calculated with the inferCNV R package of the Trinity CTAT project ("https://github.com/broadinstitute/inferCNV") and depicted as heatmap. Healthy T cells were used as reference.

To determine RNA velocities, we used the R package velocyto (13). Velocyto aims to predict future transcriptional cell states by assuming that unspliced, intronic reads can be interpreted as nascent RNAs. The ratio of unspliced and spliced reads can be used to infer RNA turnover rates (aka RNA velocities) and for estimation of future transcriptome states. These futures states of a cell can also be projected on lower dimensional embeddings to visualize transition of cell clusters. Here, transitions are visualized as arrows. Arrow origins indicate the current and arrow tips the estimated future transcriptional state. We applied the velocyto package following the R tutorial according to La Manno et al. using default settings (13). In order to quantify transition distances, we used principal component analysis (PCA) embedding. In the PCA space, arrow lengths scale linearly with the dimension and can then be associated with changes in gene expression; e.g. long arrows indicating large changes. Thus, populations with short arrows indicate stationary cells while populations with longer arrows transitioning cells. To describe transition dynamics within and across samples, we measured the transition distance (i.e., the length of the arrow) from current to future cell states in the PCA reduced representation. Transition distances were calculated by: td = (ΔPC1² + tΔPC2²)^1/2 for each cell and normalized to a value between 0 and 1 (td_normalized = (td - min(td))/max(td)). Median transition distances and the inter-quartile range (IQR: calculated as 75% quantile - 25% quantile) were used to compare dynamics among samples and variations within a sample.

After processing fresh, native tissues obtained from cutaneous MF plaques or tumor lesions of 12 patients into single-cell suspensions, 9 of these suspensions passed the predefined quality control criteria. These criteria required the suspensions to contain more than 90% viable cells, with at least 30% being CD3+ T cells. Two samples had to be excluded from further analysis as one patient also suffered from B-CLL and in the other, the initial diagnosis of CTCL was revised to a fixed toxic drug eruption ( Supplementary Figure 1A ). Histomorphologically, there is a diffuse and dense infiltration of atypical lymphocytes in the dermis. These lymphocytes have irregular, hyperchromatic, cerebriform nuclei and prominent nucleoli, often extending into the subcutaneous fat. In some instances, papillary dermal fibrosis is evident. Despite the plaque/tumor stage, most cases exhibit some degree of epidermotropism, including the formation of Pautrier’s micro-abscesses ( Supplementary Figure 1B ). Figure 1 depicts the individual disease course, including the therapies used, for each of these 7 patients. As commonly observed for cohorts of MF patients, they had widely varying disease histories, e.g., disease duration varied from less than one to 17 years and the number and type of therapies administered.

Malignant and reactive T cells freshly isolated from cutaneous MF lesions were scrutinized by combined scRNA and scTCRseq ( Figure 2A ). To ensure an unbiased analysis, we did not enrich any cells prior to single-cell analysis. Instead, we aimed to ensure that samples contained a sufficient number of malignant MF cells as well as reactive T cells. We proceeded with scRNAseq only if the single-cell suspension contained at least 90% viable cells, with at least 30% being CD3+ T cells. A lower proportion of T cells would result in an insufficient number of sequenced MF and reactive T cells, compromising the reliability of data on tumor cell heterogeneity and the T-cell microenvironment. Cells were annotated according to cell type specific marker gene expression using the Panglao Database as reference ( Supplementary Figure 2 ). From a total of 36,722 cells of all patients, 18,630 were identified as T cells based on their expression of CD3 genes. Initial analyses were focusing on T cells. Clustering of the cells and the visualization using uniform manifold approximation and projection (UMAP) revealed one large cluster mixed with cells from all patients and seven distinct clusters with cells originating from only one patient each ( Figure 2B ). Consideration of the TCR sequencing data demonstrated that the mixed large cluster consisted of polyclonal T cells, whereas the individual clusters derived from each patient consisted of only one clonal T-cell population each ( Figure 2C ). Interestingly, the cells in the distinct cluster of patient 7 expressed a clonal, but non-productive TCR ( Supplementary Figure 3 ). The proportion of the malignant MF cells was high in all patients, but varied substantially, ranging from 47% to 98% ( Figure 2D ).

Given the clear differentiation between malignant T-cell clones and the corresponding reactive T-cell infiltrate in every patient, we could then proceed to compare mRNA expression of established CTCL marker genes between these two distinct cell populations within cutaneous MF lesions ( Figure 2E ) (14). Specifically, we examined the expression of CD4, IL2RA, JAK3, STAT3, STAT5A, TOX, SELPLG, CCR4, CCR10 and CXCR3, and found that only STAT3 mRNA was expressed in levels that could be reliably detected by scRNAseq in the malignant T-cell clones of all MF samples. Nevertheless, the transcription factor TOX and the chemokine receptor CXCR3 were detected in the malignant clones of 6 of 7 patients ( Figure 2E ); genes associated with skin homing or retention, SELPLG (coding for CD162, aka CLA) and CCR10 were also expressed in most lesions, but to variable degrees. To determine whether the expression pattern of other genes may be better suited to distinguish malignant from reactive T cells by scRNAseq, we assessed differential gene expression in the pooled data from all patients ( Figure 2F ). Surprisingly, the most discriminative gene expression signature for malignant T cells was related to the mechanism or regulation of T-cell cytolytic activity such as GZMB (granzyme B), CTSW (cathepsin W, lymphopain) and GNLY (granulysin) ( Figure 2F ). Other genes with a stronger mRNA expression in MF cells were LGALS1 (galectin-1), which suppresses Th1 immune responses (15), CD9, a tetraspanin involved in T-cell activation, migration and exosome biogenesis (16), and the keratins KRT81 and KRT86, both mediating cell migration and tissue invasion ( Figure 2F ) (17). Genes with a higher mRNA expression in the benign T cells infiltrating MF lesions were members of the Krüppel-like family of transcription factors KLF2 and KLF3, as well as genes associated with T-cell memory function, including IL7R, CD27 and SELL (CD62L).

To overcome a possible bias in differential gene expression analysis by using the compiled data set of all patients (e.g., due to the variable number of malignant T cells), we next performed the analysis between reactive and malignant T cells of each patient individually. Twelve genes were differentially expressed and upregulated in the malignant cells present in four or more lesions ( Figure 2G ). Among these were pro-inflammatory genes associated with T-cell activation and differentiation, such as BATF (18), IL32 (19), MIF (20), CD9 (21), and genes involved in T cell metabolism ENO1 (enolase 1) (22) and KDSR (3-ketodihydrosphingosine reductase) (23). The AP-1 family transcription factor BATF specifically mediates the differentiation of T-helper 17 cells (Th17). IL-32, MIF and CD9 are involved in the polarization of myeloid cells, with which the malignant cells may interact via NKp30 (encoded by NCR3) to stimulate their own growth (24). Interestingly, both glycolytic enzymes ENO1 and KDSR have been shown to mediate the metabolic adaptation to hypoxic stress - a likely feature of the cutaneous microenvironment (22). The two S100 genes S100A4 and S100A6, which are involved in many cellular processes including inflammation and cancer progression, were also highly expressed in the malignant T cells (25, 26). Finally, the strong expression of mRNA encoding galectin-1 and 3 in malignant T cells suggests the secretion of these proteins; both have been shown to suppress anti-tumor responses, favor Th2 inflammation and polarize the epidermal environment to a pro-tumorigenic state (27).

Shared nearest neighbor clustering separates malignant T cells from individual patients, indicating substantial transcriptional differences between patients’ malignancies. To scrutinize differences and commonalities, we focused exclusively on the malignant T cells: the two-dimensional presentation of mRNA expression profiles by UMAP already revealed clear differences between the malignant T cells of individual patients ( Figure 3A ). Indeed, most of the genes mentioned above were expressed markedly distinct expression levels among the different patients. Figure 3B depicts a heat map of the top 10 differentially upregulated genes in MF cells between patients demonstrating a unique mRNA expression profile for each individual patient. In detail, MF cells of patient 1, 2, 5 and 6 possess a transcriptomic phenotype characteristic for memory/naïve T cells characterized by expression of IL7R, CCR7 and CD27 ( Figure 3C ). MF cells of patient 1, 2 and 5 additionally express mRNA coding for the chemoattractant CXCL13, which has been described to enhance the migratory capacity of CTCL cells (28). Patient 6’s MF cells express the cytokines IL4 and IL13, both associated with a Th2 phenotype (29). In contrast, malignant T cells of patients 3, 4 and 7 express genes associated with the regulation or execution of cytolytic activity such as GZMA, GNLY and NKG7 ( Figure 3C ); in patient 7, MF cells also strongly express GZMB and PRF1 mRNA. The high expression of glycolysis-related genes in patient 3’s tumor cells is likely caused by expression of high levels of miR-155, a microRNA promoting malignant T-cell proliferation (30). With respect to T-cell exhaustion markers, such as PDCD1 (encoding PD-1), CTLA4, HAVCR2 (encoding TIM-3), LAG3, TIGIT and CD160, which have been reported to be highly expressed in CTCL cells, we also observed a marked interpatient heterogeneity ( Figure 3D ) (31). For example, expression of TIGIT, polarizing toward Th2 immune responses, was detectable in only half of the patients (32). For LAG3, PDCD1 and CTLA4, only patient 5 had a high number of expressing cells in the cutaneous MF lesion, whereas HAVCR2 was only expressed in cells from patient 7 ( Figure 3D ).

The heterogeneity of MF cells within a patient is also reflected by the number of clusters identified by shared-nearest-neighbor clustering ( Figure 4A ) (8). However, it should be noted that the number of identified clusters also depends on the number of sequenced cells, so for samples with a low number of sequenced MF cells, e.g., patient 3 and 5, this analysis should be interpreted with caution. These patients were therefore not included in the correlation of tumor heterogeneity with RNA expression dynamics. The latter was determined using velocyto (13), which allows to infer the direction and extent of transcriptional changes from the ratio of unspliced to spliced RNA. To quantify these transcriptional changes, we calculated the median transition distance in the PC space, indicating the magnitude of transcriptional change, and the inner quartile range (IQR), indicating the variance of transcriptional changes between all cells within a sample. This analysis again revealed large differences in the dynamics of mRNA processing between patients, thus providing additional evidence toward the heterogeneity within MF cells ( Figure 4B ). For depiction of cell transitions, we projected RNA velocities on the UMAP and used vector fields to combine individual cell velocities into neighboring group velocities ( Figure 4C ). These visualizations demonstrate that, for example, the highly diverse malignant T-cell population in patient 4 appears to be transcriptionally rather static, whereas the less diverse MF cell population in patient 2 shows considerable transcriptional dynamics.

The observed transcriptomic differences between and within individual patients may be due to the previously reported genetic heterogeneity observed in late-stage MF; indeed, only few common driver genes were identified (3, 8, 33). To test this notion, we derived copy number variations (CNVs) from the single-cell transcriptome data using the inferCNV algorithm. Sequencing data from the benign reactive T cells served as the normal reference. This analysis demonstrated that each malignant T-cell clone possessed an individual CNV pattern, which confirms the previous reports ( Figure 5A ). Of note, we observed more genetic changes in patients whose long-standing MF was treated with multiple therapies, e.g., patients 5 and 6; this suggests that the duration of the disease and the number of therapies, especially ionizing radiation, may result in more genetic changes ( Figure 1 ). Amplifications were frequently found on chromosome 1p (ENO1, TNFRSF1B; n=5); deletions on chromosome 6 (MHCI and MHCII, n=4). However, the characteristic deletions of chromosomes 10q (PTEN, FAS, NFKB2) and 17p (TP53) reported for leukemic MF and SS were only observed in patient 5 and 2, respectively. The observed transcriptional heterogeneity of MF cells within individual patients was also reflected to some extent in the inferred CNV signature, as shown by the example of patients 2 and 5 ( Figures 5B, C ).