Identification of a Distinct miRNA Regulatory Network in the Tumor Microenvironment of Transformed Mycosis Fungoides

Simple Summary Transformed mycosis fungoides (LCT-MF) is a histopathological marker of poor prognosis and associated with worse survival. We compared miRNA and mRNA expression profiles of LCT-MF with classic MF and found a distinct miRNA regulatory network modulated immunosuppressive tumor microenvironment in LCT-MF. Our findings provide novel insights and therapeutic targets for LCT-MF. Abstract Large cell transformation of mycosis fungoides (LCT-MF) occurs in 20–50% of advanced MF and is generally associated with poor response and dismal prognosis. Although different mechanisms have been proposed to explain the pathogenesis, little is known about the role of microRNAs (miRs) in transcriptional regulation of LCT-MF. Here, we investigated the miR and mRNA expression profile in lesional skin samples of patients with LCT-MF and non-LCT MF using RNA-seq analysis. We found miR-146a and miR-21 to be significantly upregulated, and miR-708 the most significantly downregulated miR in LCT-MF. Integration of miR and mRNA expression profiles revealed the miR-regulated networks in LCT-MF. Ingenuity pathway analysis (IPA) demonstrated the involvement of genes for ICOS-ICOSL, PD1-PDL1, NF-κB, E2F transcription, and molecular mechanisms of cancer signaling pathways. Quantitative real time (qRT)-PCR results of target genes were consistent with the RNA-seq data. We further identified the immunosuppressive tumor microenvironment (TME) in LCT-MF. Moreover, our data indicated that miR-146a, -21 and -708 are associated with the immunosuppressive TME in LCT-MF. Collectively, our results suggest that the key LCT-MF associated miRs and their regulated networks may provide insights into its pathogenesis and identify promising targets for novel therapeutic strategies.


Introduction
Cutaneous T cell lymphoma (CTCL) is a disfiguring, incurable malignancy with a poor prognosis for those with advanced stage disease. Mycosis fungoides (MF) is the most common primary CTCL and is characterized by expansion of malignant T-cells in a chronic inflammatory environment in the skin [1]. MF evolves through patches and plaques to tumors and/or erythroderma with or without systemic involvement. Patients with early-stage disease have a 10-year survival ranging between 64-98%, whereas patients with advanced stage disease have a less favorable prognosis, with a 10-year survival of 20-42% [2,3]. Large cell transformation (LCT) of CTCL, histopathologically defined by sheets or microscopic nodules of atypical T cells with a large cell phenotype four times larger than normal lymphocytes that make up at least 25% of the dermal infiltrate, is associated with rapidly progressing disease and dismal prognosis and confers a median survival of 1-4 years after transformation, which is not overcome by any regimen to date [4][5][6].
Several efforts have been made to identify biomarkers of transformation that could help prompt an early and targeted treatment to lead to a more favorable clinical outcome. The most promising targets of these investigations are the genetic aberrations and epigenetic changes that have been associated with disease progression and survival in patients with MF [7]. Notably, many genetic and chromosomal abnormalities have been identified in the pathogenesis and progression of MF, such as amplification of JUNB, which is involved in T-cell proliferation, differentiation, and apoptosis, [8,9] and activation of the STAT family transcription factors [10]. Furthermore, decreased expression of regulating proteins p14, p15, and p16, which can interact with cyclin-dependent kinases and induce cell cycle arrest, has been described in MF [11]. Nevertheless, these genetic and chromosomal abnormalities are not sufficient to predict the outcome of LCT-MF and additional analyses are critically needed to better dissect the molecular basis of transformation.
MicroRNAs (miRs) are a class of short noncoding RNAs that regulate gene expression. They are involved in physiological processes and disease development and have an important role in DNA damage response [12]. Dysregulated miRNAs have been linked to cancer, acting either as oncogenes or tumor suppressors, and appear to be involved in the pathogenesis and progression of CTCL [13,14]. In MF, several studies have identified upregulated miR-155 in cutaneous lesions of MF, which has led to the clinical development of a miR-155 inhibitor in CTCL and other hematological malignancies [15]. The upregulation of miR-155 with concomitant downregulation of miR-203 and miR-205 in MF was found to be useful to discriminate MF from benign dermatoses [16]. Notably, Lindahl et al. [17] developed a 3-miR classifier, based on miR-106b-5p, miR-148a-3p, and miR-338-3p, that appears to predict disease progression of early MF, and Manso et al. demonstrated that miRs regulate MF progression [18]. MiR expression profiling has also been used for prognostication in tumor stage and subtypes of MF; high expression of miR-155 and low expression of miR-200b correlated with worse prognosis [13]. However, the mechanisms by which the miRs contribute to the signaling pathways and microenvironmental immune cell profiles in LCT-MF remain unknown. The goal of this study was to identify differentially expressed miRs, the regulatory network targeting pathways of LCT-MF, and their effects on the TME, providing insights into pathogenesis and potentially identifying promising targets for novel therapeutic strategies in LCT-MF.

Patients and Tissue Samples
A prospective consecutive case-series of 28 patients (20 males, 8 females; median age: 63.6 years, range: 30.8-85.4 years) diagnosed with MF received little to no treatment prior to their initial visit in our multidisciplinary cutaneous lymphoma clinic. Subjects had signed informed consent for specimen collection for research purposes, in accordance with the Declaration of Helsinki and approval by the City of Hope Institutional Review Board. Patients were diagnosed according to the revised 2018 WHO-EORTC classification [2] and classified according to the revised staging system for MF [19]. Table 1 summarizes patient characteristics. Skin disease was classified according to T stage. T1 and T2 are defined by plaques involving <10% and ≥10% of the body surface area, respectively, whereas T3 is defined by ≥1 tumor and T4 by erythroderma. Relevant demographic and clinical information was collected from the electronic medical records including age, gender, race/ethnicity, diagnosis, stage, histopathologic findings such as transformation status, and assessment of the skin tumor burden by the severity-weighted assessment tool (mSWAT). LCT-MF is defined as showing large cells (4 times the size of a small lymphocyte) in 25% or more of the dermal infiltrate or forming microscopic nodules [6].

RNA Sequencing and Library Preparation
Total RNA from 34 paraffin-embedded (FFPE) lesional CTCL skin biopsy sections (size uniform, plaques and tumors ≥ 1 cm 2 ), which included 14 samples of large cell transformation (7 plaques, 7 tumors) and 20 non-transformed lesions (12 plaques, 8 tumors) was extracted with a miRNAeasy FFPE kit (Qiagen, Valencia, CA, USA). Library preparation and RNA sequencing were performed as described in our previous publication [20]. Briefly, 300 ng of the total RNA isolation was used. Small RNA library preparation was performed following the Illumina protocol with minor modifications. The library was created from each RNA sample by 3 adapter ligation, 5 RT primer annealing, 5 adapter ligation, reverse transcription, and PCR amplification. Sequencing of 50 cycles was performed on a HiSeq 2500 (Illumina, San Diego, CA, USA). Demultiplexing of the raw sequencing data and generation of the fastq files were performed using the Real-Time Analysis (RTA) software. For Ribo-Zero RNA Sequencing, total RNA was rRNA-depleted using the Ribo-Zero low input kit for Human/Mouse/Rat (Kapa Biosystems, Wilmington, NC, USA). RNA-seq libraries were prepared from Ribo-Zero mRNA-enriched material using KAPA Stranded RNA-Seq Library Preparation Kit (Illumina Platforms) (Kapa Biosystems, Wilmington, NC, USA), using 10 cycles of PCR amplification. Libraries were purified using AxyPrep Mag PCR Clean-up kit (Thermo Fisher Scientific, Irwindale, CA, USA). Sequencing was performed on an Illumina ® Hiseq 2500 (Illumina, San Diego, CA, USA) instrument using the TruSeq SR Cluster Kit V4-cBot-HS (Illumina ® ) to generate 51 bp single-end reads sequencing with v4 chemistry. Quality control of RNA-Seq reads was performed using FastQC.

Sequencing Data Analysis
miR and mRNA-seq data were analyzed as previously described [20] with the following modification: The mapping table was created using the human miR mature sequences from miRBase release18 and aligned back to human hg19 genome afterwards. Differential expression analysis was conducted using Bioconductor package "edgeR", with p value < 0.01 and fold change > 1.5 as the cutoff values for miRNA data. For RNA-seq data, "limma" was used on log2 transformed RPKM values, with p value < 0.05 and fold change > 1.5 as the cutoff values. Hallmark pathway enrichment analyses were carried out using gene set enrichment analysis (GSEA) with genes sorted by-log10 (p value) with sign determined by the fold change between MF-LCT and non-LCT samples. An FDR ≤ 0.05 and the involvement of at least 2 genes were used as the cutoff criteria. These data sets are deposited in the Gene Expression Omnibus (GEO) repository under the accession numbers GSE181130 and GSE113113.

Functional Gene Set Enrichment Analysis
An interaction network of mRNAs in LCT-MF that are upregulated and associated with miR-146a and -21 was generated using Ingenuity Pathway Analysis (IPA, Ingenuity ® Systems, www.ingenuity.com, accessed on 10 October 2020) with fold change > 1.5, FDR < 0.05.

Statistical Analysis
Kaplan-Meier survival curves with Mantel-Cox test was used for survival analysis. Student's t test was used to analyze the data when not otherwise specified. All p values were two-sided, and the statistical significance level was set at p < 0.05. Correlation between differentially expressed miRs and immune cell gene score was analyzed by Pearson's correlation coefficient using Bioconductor package "edgeR". Statistical analyses were carried out using Prism 8.2.0 (GraphPad, La Jolla, CA, USA) when not otherwise specified.

Outcomes for Patients with LCT-MF Are Discouraging
Clinical stage is an important determinant of risk for disease progression (RDP) and overall survival (OS) [3]. The clinical characteristics and distribution of patients are summarized in Table 1 Figure 1C. On histopathology, an atypical epidermotropic lymphoid infiltrate with large cell morphology is noted. The large pleomorphic cells uniformly show high expression for PD1 and PD-L1 ( Figure 1D).

The miR Expression Profile in LCT-MF Is Distinct from That of Non-LCT
We determined the miR profiles from skin biopsies of 19 LCT-MF and 15 non-LCT patients using small RNA sequencing. The unsupervised heatmap of miR analysis reveals different miR clusters for LCT and non-LCT-MF (Figure 2A). Furthermore, comparison of the miR expression levels in LCT-MF plaques and tumors with that in non-LCT-MF identified 17 differentially expressed miRs with an adjusted p value < 0.05 ( Table 2). The highest-ranked miRNAs included miR-21-3p (abbreviated as miR-21 below), miR-146a-3p (abbreviated as miR-146a below), miR-136-5p, miR-889 and miR-539-3p. The majority (seven miRs) were downregulated, with a log fold change (log FC) < −2. The top downregulated miRs included miR-708-5p, miR-5701, and miR-3653, which are widely viewed as tumor suppressor miRs and noted in many cancers during tumor progression. The upregulation of miR-21 and miR-146a and downregulation of miR-708 in LCT-MF were verified by qRT-PCR analysis. As shown in Figure 2B, both miR-146a (p < 0.05) and -21 (p < 0.05) expression levels were found to be elevated and miR-708 (p < 0.05) was decreased in LCT-MF, validating the findings of RNA-seq analysis. Thus, these results suggested that miR-146a and -21 might be critical for LCT-MF progression.  The 2ˆ(-delta delta CT) method was used as a relative quantification strategy for data analysis. Results are shown as means ± SD, and differences were tested for significance using Student's t test (* p < 0.05).

MRNA Expression Profiling Distinguishes LCT-MF from Non-LCT MF
We compared the global mRNA expression profile of LCT vs. non-LCT cases chosen by histopathologic criteria. Hierarchical clustering revealed two distinct clusters segregated by transformation status, shown in the unsupervised mRNA heatmap ( Figure 3A). The left cluster reveals LCT-MF cases (plaques and tumors); the right cluster identifies non-LCT MF. To identify the potential functions of the identified differentially expressed genes (DEGs), we performed Hallmark pathway analysis. The top enriched Hallmark pathway terms for both miRs included E2F targets, MYC targets, G2M checkpoint, oxidative phosphorylation, DNA repair and MTORC1. The minimally enriched and negatively enriched Hallmark pathway terms among DEGs are shown in Figure 3B. qRT-PCR was performed to verify the expression patterns of mRNAs determined by RNA seq. As shown in Figure 3C, the findings showed that the expression patterns of STAT3 (p < 0.05), NF-κB p65 (p < 0.05) and TGFβ1 (p < 0.05) were consistent with the RNA-seq data. These findings indicated that in most cases the qRT-PCR results were consistent with those of RNA-seq, which implied that the mRNA expression profile in LCT-MF is reliable.

Distinct Microenvironmental Immune Cell Populations in LCT-MF
To understand the diversity of immune cell types that drive LCT-MF pathogenesis, we profiled microenvironmental immune cell gene signatures with CIBERSORT. The analysis of immune cell abundance estimated using CIBERSORT identified signatures for naïve CD4+ T cells (p < 0.01), resting dendritic cells (p < 0.001), gamma delta T cells (p < 0.05), M0 macrophages (p < 0.05), M2 macrophages (p < 0.05) as more abundant, whereas CD4+ T cells (p < 0.001) and mast cells were less abundant in LCT-MF than in non-LCT ( Figure S1A). When comparing the PD1 and PD-L1 transcripts in LCT-MF with those in non-LCT, we found a higher frequency of PD1 transcripts in CD8+ T cells and CD4+ T cells in LCT-MF than in non-LCT-MF and PD1 transcripts in neutrophils in non-LCT over that in LCT; a higher frequency of PD-L1 transcripts were detected in B cells, CD8+ T cells, and eosinophils in non-LCT than in LCT, whereas neutrophil expression was relatively higher in LCT ( Figure S1B). Furthermore, sequencing analysis was performed using gene sets representing B cells, CD8+ tumor-infiltrating lymphocytes (TILs), tumor-associated macrophages (TAMs), dendritic cells (DC), myeloid-derived suppressor cells (MDSC), regulatory T cells (Tregs), and Th1/Th2 cell memory phenotypes to analyze the immune cell gene profile in the CTCL microenvironment of LCT-MF versus non-LCT. Heatmaps were generated using the defined genes representing immune cell types ( Figure 5A,B). We used the Kruskal-Wallis test to assess the average expression scores of relevant immune cell related genes including TAMs (M2), MDSCs, iTregs, B cells, and mature DCs stratified by transformation status ( Figure 5C). Our data indicate that microenvironmental immune cell gene profiling in LCT-MF demonstrated a significant elevated gene signature for inducible (i) Tregs (p = 0.036), PD-L1+ M2 macrophages (p = 0.016), and PD1+ M2-macrophages (p = 0.013) and a trend for elevated genes for MDSCs (p = 0.067) and for a decreased gene signature for mature DCs (p = 0.1), which may indicate an immune cell profile associated with immunosuppression and progression of LCT. There was no significant difference seen in B cell -related genes (p = 0.82) between LCT-MF vs. MF non-LCT. The findings indicate an immunosuppressive TME in LCT-MF.

MiR Expression Correlates with the Immune Cell Gene Score in LCT-MF
Lastly, we wanted to investigate the relationship between miR expression and the immune cell gene signature in LCT-MF. To this end, a Pearson's correlation coefficient analysis was carried out between miR expression and the average gene score of immune cells. miR-21, which is highly expressed in LCT-MF, showed a strong positive correlation with iTreg (r = 0.46, p = 0.0073) and MDSC (r = 0.37, p = 0.037) gene expression ( Figure 6A); the upregulated miR-146a also positively correlated with iTreg gene score (r = 0.6, p = 0.00025) in LCT-MF ( Figure 6B). Although not significant, there was a trend toward a positive correlation with MDSCs. Moreover, the downregulated miR-708 negatively correlated with the TAM (M2) gene score (r = 0.46, p = 0.0065; Figure 6C), suggesting miR-708-mediated effects could inhibit an immunosuppressive environment, through inhibition of TAM (M2) gene expression. A trend was also seen for negative correlation with iTregs and MDSCs and a positive correlation with mature DCs. In summary, upregulated miR-21 and miR-146a and downregulated miR-708 drive an immunosuppressive tumor microenvironment that may contribute to the transformed phenotype in MF and impact tumor growth.

Discussion
The role of epigenetic dysregulation in the pathogenesis of LCT-MF is an area of active research and clinical interest. Notably, a phase 1 study of MRG-106 (cobomarsen), an inhibitor of miR-155, has shown efficacy in MF patients. The role and regulation of miR profiles in LCT-MF pathogenesis has not previously been explored. Therefore, the identification of the miRs targeting regulatory signaling pathways and their effects on the immune microenvironment in LCT-MF is critical to develop effective treatment strategies.
Previous studies have identified LCT as an independent prognostic factor in patients with MF, portending a poor prognosis and outcome [21]. Our findings confirm previously published data [21,22] highlighting a dismal 2-year OS of 14.3% in patients with LCT-MF versus 80% for non-LCT MF. Transformation occurs predominantly in advanced/progressing stages of MF presenting with cutaneous and/or extracutaneous tumors. However, in our cohort, a subset of early stage LCT-MF, clinically presenting with cutaneous plaque, was included. We found 12 miRNAs (five upregulated and seven downregulated) differentially expressed in LCT-MF tumors and plaques compared to non-LCT MF. The signatures in plaques were similar to tumors and may indicate early prognostic factors for transformation and progression. In non-LCT-MF, we had recently established that the positive correlation of miR-155, -130, and-21 with an immune checkpoint (IC) expression profile provides justification for miRs as a putative therapeutic target to reverse T cell exhaustion in CTCL [23]. To our surprise, the upregulated miRs in LCT-MF cases did not include miR-155. Although data show that miR-155 is linked to the oncogenesis of MF, with higher expression noted in tumor lesions [24], our findings suggest that miR-155 regulation is not critical for transformation. Among the upregulated miRs in LCT-MF, we found that the miR with the highest fold difference was miR-146a. MiR-146a has emerged as key regulator of immune responses and tumorigenesis in various cancers [25,26]. MiR-21 was previously identified as the top related miR associated with disease progression and as highly upregulated in advanced CTCL stages [27]. Lindahl LM et al. found miR-21 expression was driven by JAK3/STAT5 in both malignant and non-malignant T cells in MF [28].
To date, miR-regulated signaling pathways remain largely unknown in LCT-MF. In this study, we identified the upregulated miR-related signaling network in LCT-MF including ICOS-ICOSL signaling, PD1-PDL1 signaling, NF-κB signaling, E2F transcription, and molecular mechanisms of cancer pathways. Although these pathways have been consistent with previous studies in MF, other specific pathways of LCT-MF in this analysis, such as E2F transcription, have not been previously reported and are worthwhile to explore, as the interaction of malignant cells with E2F, for one, promotes a microenvironment suitable for tumor cell growth, survival, and aggressiveness [29,30]. Among the downregulated miRs, the miR with the most statistically significant log FC was the hsa-miR-708. There are opposing results among studies as to the oncogenic or tumor suppressor function of miR-708, which may be due to the identity of targets driving various cancers [31]. Little is known in T cell malignancies. Low miR-708 expression was found to correlate with high-risk T cell progenitor-derived T-ALL and poor prognosis [32]. Notably, miR-708 may function as a tumor suppressor by directly targeting the innate checkpoint CD47 in T-ALL [33]. Restoration of miR-708 expression in a T-ALL cell line promoted phagocytosis by macrophages via downregulation of CD47 to eliminate T-ALL cells. Considering the emerging role of miRs involved in immune responses, these miR candidates seem to be of great interest for further exploration in LCT-MF.
The importance of the cellular microenvironment in various solid and hematologic malignancies including CTCL has been increasingly recognized [3,34,35]. We sought to identify the effects of miRs on the immune cell gene profile in MF-LCT. A significant gene score was noted for iTregs and PD-L1+ M2 in MF-LCT when compared to non-LCT, suggesting that iTregs and PD-L1+ M2 maybe the most abundant cell types in the tumor microenvironment in localized LCT-MF. Although not statistically significant, possibly due to the overall small sample size, there was a trend toward statistically significant increased TAM (M2) MDSC and B cell gene scores in LCT-MF and a decreased mature DC gene score. Tumor-derived miRs may mediate intercellular communication affecting the maturation status of DCs. Several studies demonstrated a decrease in DC subsets and impaired cell-mediated immunity in MF [36,37]. Manso et al. found that miR-146a was highly expressed in advanced stages of MF and directly targeted the FOXP3 gene, involved in Treg differentiation and MF progression [18]. We and others have shown the pro-tumorigenic role of TAM2-macrophages in MF [38]. In line with these observations, our findings support an immunosuppressive TME in LCT-MF. Recent studies have shown that miRs play a role in remodeling the TME [39,40]; however, it remains unknown which miRs regulate the immunosuppressive microenvironment in LCT-MF. We identified specific miR signatures associated with immune cell mRNA signatures. Notably, upregulated miR-21 and miR-146a positively correlated with iTreg gene score, and downregulated miR-708 negatively corelated with TAM (M2) gene score in MF-LCT. These findings show that not only the upregulated miRs, but also the downregulated miRs, promote the formation of an immunosuppressive microenvironment in LCT-MF.

Conclusions
In conclusion, our study identified miR-regulated signaling networks and their effects on immune cell gene profiles in LCT-MF, which may prove relevant to pathogenesis and hold promise for the identification of potential therapeutic targets in this disease.