SCENIC+: single-cell multiomic inference of enhancers and gene regulatory networks.

Joint profiling of chromatin accessibility and gene expression in individual cells provides an opportunity to decipher enhancer-driven gene regulatory networks (GRNs). Here we present a method for the inference of enhancer-driven GRNs, called SCENIC+. SCENIC+ predicts genomic enhancers along with candidate upstream transcription factors (TFs) and links these enhancers to candidate target genes. To improve both recall and precision of TF identification, we curated and clustered a motif collection with more than 30,000 motifs. We benchmarked SCENIC+ on diverse datasets from different species, including human peripheral blood mononuclear cells, ENCODE cell lines, melanoma cell states and Drosophila retinal development. Next, we exploit SCENIC+ predictions to study conserved TFs, enhancers and GRNs between human and mouse cell types in the cerebral cortex. Finally, we use SCENIC+ to study the dynamics of gene regulation along differentiation trajectories and the effect of TF perturbations on cell state. SCENIC+ is available at scenicplus.readthedocs.io .

A cellular hierarchy in melanoma uncouples growth and metastasis.

Although melanoma is notorious for its high degree of heterogeneity and plasticity1,2, the origin and magnitude of cell-state diversity remains poorly understood. Equally, it is unclear whether growth and metastatic dissemination are supported by overlapping or distinct melanoma subpopulations. Here, by combining mouse genetics, single-cell and spatial transcriptomics, lineage tracing and quantitative modelling, we provide evidence of a hierarchical model of tumour growth that mirrors the cellular and molecular logic underlying the cell-fate specification and differentiation of the embryonic neural crest. We show that tumorigenic competence is associated with a spatially localized perivascular niche, a phenotype acquired through an intercellular communication pathway established by endothelial cells. Consistent with a model in which only a fraction of cells are fated to fuel growth, temporal single-cell tracing of a population of melanoma cells with a mesenchymal-like state revealed that these cells do not contribute to primary tumour growth but, instead, constitute a pool of metastatic initiating cells that switch cell identity while disseminating to secondary organs. Our data provide a spatially and temporally resolved map of the diversity and trajectories of melanoma cell states and suggest that the ability to support growth and metastasis are limited to distinct pools of cells. The observation that these phenotypic competencies can be dynamically acquired after exposure to specific niche signals warrant the development of therapeutic strategies that interfere with the cancer cell reprogramming activity of such microenvironmental cues.

Tyrosine-Dependent Phenotype Switching Occurs Early in Many Primary Melanoma Cultures Limiting Their Translational Value.

The use of patient-derived primary cell cultures in cancer preclinical assays, including drug screens and genotoxic studies, has increased in recent years. However, their translational value is constrained by several limitations, including variability that can be caused by the culture conditions. Here, we show that the medium composition commonly used to propagate primary melanoma cultures has limited their representability of their tumor of origin and their cellular plasticity, and modified their sensitivity to therapy. Indeed, we established and compared cultures from different melanoma patients propagated in parallel in low-tyrosine (Ham's F10) or in high-tyrosine (Ham's F10 supplemented with tyrosine or RPMI1640 or DMEM) media. Tyrosine is the precursor of melanin biosynthesis, a process particularly active in differentiated melanocytes and melanoma cells. Unexpectedly, we found that the high tyrosine concentrations promoted an early phenotypic drift towards either a mesenchymal-like or senescence-like phenotype, and prevented the establishment of cultures of melanoma cells harboring differentiated features, which we show are frequently present in human clinical biopsies. Moreover, the invasive phenotype emerging in these culture conditions appeared irreversible and, as expected, associated with intrinsic resistance to MAPKi. In sharp contrast, differentiated melanoma cell cultures retained their phenotypes upon propagation in low-tyrosine medium, and importantly their phenotypic plasticity, a key hallmark of melanoma cells. Altogether, our findings underline the importance of culturing melanoma cells in low-tyrosine-containing medium in order to preserve their phenotypic identity of origin and cellular plasticity.

Analysis of long and short enhancers in melanoma cell states.

Understanding how enhancers drive cell-type specificity and efficiently identifying them is essential for the development of innovative therapeutic strategies. In melanoma, the melanocytic (MEL) and the mesenchymal-like (MES) states present themselves with different responses to therapy, making the identification of specific enhancers highly relevant. Using massively parallel reporter assays (MPRAs) in a panel of patient-derived melanoma lines (MM lines), we set to identify and decipher melanoma enhancers by first focusing on regions with state-specific H3K27 acetylation close to differentially expressed genes. An in-depth evaluation of those regions was then pursued by investigating the activity of overlapping ATAC-seq peaks along with a full tiling of the acetylated regions with 190 bp sequences. Activity was observed in more than 60% of the selected regions, and we were able to precisely locate the active enhancers within ATAC-seq peaks. Comparison of sequence content with activity, using the deep learning model DeepMEL2, revealed that AP-1 alone is responsible for the MES enhancer activity. In contrast, SOX10 and MITF both influence MEL enhancer function with SOX10 being required to achieve high levels of activity. Overall, our MPRAs shed light on the relationship between long and short sequences in terms of their sequence content, enhancer activity, and specificity across melanoma cell states.

BNIP3 promotes HIF-1α-driven melanoma growth by curbing intracellular iron homeostasis.

BNIP3 is a mitophagy receptor with context-dependent roles in cancer, but whether and how it modulates melanoma growth in vivo remains unknown. Here, we found that elevated BNIP3 levels correlated with poorer melanoma patient's survival and depletion of BNIP3 in B16-F10 melanoma cells compromised tumor growth in vivo. BNIP3 depletion halted mitophagy and enforced a PHD2-mediated downregulation of HIF-1α and its glycolytic program both in vitro and in vivo. Mechanistically, we found that BNIP3-deprived melanoma cells displayed increased intracellular iron levels caused by heightened NCOA4-mediated ferritinophagy, which fostered PHD2-mediated HIF-1α destabilization. These effects were not phenocopied by ATG5 or NIX silencing. Restoring HIF-1α levels in BNIP3-depleted melanoma cells rescued their metabolic phenotype and tumor growth in vivo, but did not affect NCOA4 turnover, underscoring that these BNIP3 effects are not secondary to HIF-1α. These results unravel an unexpected role of BNIP3 as upstream regulator of the pro-tumorigenic HIF-1α glycolytic program in melanoma cells.

Interpretation of allele-specific chromatin accessibility using cell state-aware deep learning.

Genomic sequence variation within enhancers and promoters can have a significant impact on the cellular state and phenotype. However, sifting through the millions of candidate variants in a personal genome or a cancer genome, to identify those that impact cis-regulatory function, remains a major challenge. Interpretation of noncoding genome variation benefits from explainable artificial intelligence to predict and interpret the impact of a mutation on gene regulation. Here we generate phased whole genomes with matched chromatin accessibility, histone modifications, and gene expression for 10 melanoma cell lines. We find that training a specialized deep learning model, called DeepMEL2, on melanoma chromatin accessibility data can capture the various regulatory programs of the melanocytic and mesenchymal-like melanoma cell states. This model outperforms motif-based variant scoring, as well as more generic deep learning models. We detect hundreds to thousands of allele-specific chromatin accessibility variants (ASCAVs) in each melanoma genome, of which 15%-20% can be explained by gains or losses of transcription factor binding sites. A considerable fraction of ASCAVs are caused by changes in AP-1 binding, as confirmed by matched ChIP-seq data to identify allele-specific binding of JUN and FOSL1. Finally, by augmenting the DeepMEL2 model with ChIP-seq data for GABPA, the TERT promoter mutation, as well as additional ETS motif gains, can be identified with high confidence. In conclusion, we present a new integrative genomics approach and a deep learning model to identify and interpret functional enhancer mutations with allelic imbalance of chromatin accessibility and gene expression.

A Multi-Omics Analysis of Metastatic Melanoma Identifies a Germinal Center-Like Tumor Microenvironment in HLA-DR-Positive Tumor Areas.

The emergence of immune checkpoint inhibitors has dramatically changed the therapeutic landscape for patients with advanced melanoma. However, relatively low response rates and a high incidence of severe immune-related adverse events have prompted the search for predictive biomarkers. A positive predictive value has been attributed to the aberrant expression of Human Leukocyte Antigen-DR (HLA-DR) by melanoma cells, but it remains unknown why this is the case. In this study, we have examined the microenvironment of HLA-DR positive metastatic melanoma samples using a multi-omics approach. First, using spatial, single-cell mapping by multiplexed immunohistochemistry, we found that the microenvironment of HLA-DR positive melanoma regions was enriched by professional antigen presenting cells, including classical dendritic cells and macrophages, while a more general cytotoxic T cell exhaustion phenotype was present in these regions. In parallel, transcriptomic analysis on micro dissected tissue from HLA-DR positive and HLA-DR negative areas showed increased IFNγ signaling, enhanced leukocyte adhesion and mononuclear cell proliferation in HLA-DR positive areas. Finally, multiplexed cytokine profiling identified an increased expression of germinal center cytokines CXCL12, CXCL13 and CCL19 in HLA-DR positive metastatic lesions, which, together with IFNγ and IL4 could serve as biomarkers to discriminate tumor samples containing HLA-DR overexpressing tumor cells from HLA-DR negative samples. Overall, this suggests that HLA-DR positive areas in melanoma attract the anti-tumor immune cell infiltration by creating a dystrophic germinal center-like microenvironment where an enhanced antigen presentation leads to an exhausted microenvironment, nevertheless representing a fertile ground for a better efficacy of anti-PD-1 inhibitors due to simultaneous higher levels of PD-1 in the immune cells and PD-L1 in the HLA-DR positive melanoma cells.

Cross-species analysis of enhancer logic using deep learning.

Deciphering the genomic regulatory code of enhancers is a key challenge in biology because this code underlies cellular identity. A better understanding of how enhancers work will improve the interpretation of noncoding genome variation and empower the generation of cell type-specific drivers for gene therapy. Here, we explore the combination of deep learning and cross-species chromatin accessibility profiling to build explainable enhancer models. We apply this strategy to decipher the enhancer code in melanoma, a relevant case study owing to the presence of distinct melanoma cell states. We trained and validated a deep learning model, called DeepMEL, using chromatin accessibility data of 26 melanoma samples across six different species. We show the accuracy of DeepMEL predictions on the CAGI5 challenge, where it significantly outperforms existing models on the melanoma enhancer of IRF4 Next, we exploit DeepMEL to analyze enhancer architectures and identify accurate transcription factor binding sites for the core regulatory complexes in the two different melanoma states, with distinct roles for each transcription factor, in terms of nucleosome displacement or enhancer activation. Finally, DeepMEL identifies orthologous enhancers across distantly related species, where sequence alignment fails, and the model highlights specific nucleotide substitutions that underlie enhancer turnover. DeepMEL can be used from the Kipoi database to predict and optimize candidate enhancers and to prioritize enhancer mutations. In addition, our computational strategy can be applied to other cancer or normal cell types.

Robust gene expression programs underlie recurrent cell states and phenotype switching in melanoma.

Melanoma cells can switch between a melanocytic and a mesenchymal-like state. Scattered evidence indicates that additional intermediate state(s) may exist. Here, to search for such states and decipher their underlying gene regulatory network (GRN), we studied 10 melanoma cultures using single-cell RNA sequencing (RNA-seq) as well as 26 additional cultures using bulk RNA-seq. Although each culture exhibited a unique transcriptome, we identified shared GRNs that underlie the extreme melanocytic and mesenchymal states and the intermediate state. This intermediate state is corroborated by a distinct chromatin landscape and is governed by the transcription factors SOX6, NFATC2, EGR3, ELF1 and ETV4. Single-cell migration assays confirmed the intermediate migratory phenotype of this state. Using time-series sampling of single cells after knockdown of SOX10, we unravelled the sequential and recurrent arrangement of GRNs during phenotype switching. Taken together, these analyses indicate that an intermediate state exists and is driven by a distinct and stable 'mixed' GRN rather than being a symbiotic heterogeneous mix of cells.

The EMT Transcription Factor ZEB2 Promotes Proliferation of Primary and Metastatic Melanoma While Suppressing an Invasive, Mesenchymal-Like Phenotype.

Epithelial-to-mesenchymal transition (EMT)-inducing transcription factors (TF) are well known for their ability to induce mesenchymal states associated with increased migratory and invasive properties. Unexpectedly, nuclear expression of the EMT-TF ZEB2 in human primary melanoma has been shown to correlate with reduced invasion. We report here that ZEB2 is required for outgrowth for primary melanomas and metastases at secondary sites. Ablation of Zeb2 hampered outgrowth of primary melanomas in vivo, whereas ectopic expression enhanced proliferation and growth at both primary and secondary sites. Gain of Zeb2 expression in pulmonary-residing melanoma cells promoted the development of macroscopic lesions. In vivo fate mapping made clear that melanoma cells undergo a conversion in state where ZEB2 expression is replaced by ZEB1 expression associated with gain of an invasive phenotype. These findings suggest that reversible switching of the ZEB2/ZEB1 ratio enhances melanoma metastatic dissemination. SIGNIFICANCE: ZEB2 function exerts opposing behaviors in melanoma by promoting proliferation and expansion and conversely inhibiting invasiveness, which could be of future clinical relevance. GRAPHICAL ABSTRACT: http://cancerres.aacrjournals.org/content/canres/80/14/2983/F1.large.jpg.

Low expression of pro-apoptotic proteins Bax, Bak and Smac indicates prolonged progression-free survival in chemotherapy-treated metastatic melanoma.

Despite the introduction of novel targeted therapies, chemotherapy still remains the primary treatment for metastatic melanoma in poorly funded healthcare environments or in case of disease relapse, with no reliable molecular markers for progression-free survival (PFS) available. As chemotherapy primarily eliminates cancer cells by apoptosis, we here evaluated if the expression of key apoptosis regulators (Bax, Bak, Bcl-2, Bcl-xL, Smac, Procaspase-9, Apaf-1, Procaspase-3 and XIAP) allows prognosticating PFS in stage III/IV melanoma patients. Following antibody validation, marker expression was determined by automated and manual scoring of immunohistochemically stained tissue microarrays (TMAs) constructed from treatment-naive metastatic melanoma biopsies. Interestingly and counter-intuitively, low expression of the pro-apoptotic proteins Bax, Bak and Smac indicated better prognosis (log-rank p < 0.0001, p = 0.0301 and p = 0.0227 for automated and p = 0.0422, p = 0.0410 and p = 0.0073 for manual scoring). These findings were independently validated in the cancer genome atlas (TCGA) metastatic melanoma cohort (TCGA-SKCM) at transcript level (log-rank p = 0.0004, p = 0.0104 and p = 0.0377). Taking expression heterogeneity between the markers in individual tumour samples into account allowed defining combinatorial Bax, Bak, Smac signatures that were associated with significantly increased PFS (p = 0.0002 and p = 0.0028 at protein and transcript level, respectively). Furthermore, combined low expression of Bax, Bak and Smac allowed predicting prolonged PFS (> 12 months) on a case-by-case basis (area under the receiver operating characteristic curve (ROC AUC) = 0.79). Taken together, our results therefore suggest that Bax, Bak and Smac jointly define a signature with potential clinical utility in chemotherapy-treated metastatic melanoma.

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

In melanoma, the lymphocytic infiltrate is a prognostic parameter classified morphologically into 'brisk', 'non-brisk' and 'absent' entailing a functional association that has never been proved. Recently, it has been shown that lymphocytic populations can be very heterogeneous, and that anti-PD-1 immunotherapy supports activated T cells. Here, we characterize the immune landscape in primary melanoma by high-dimensional single-cell multiplex analysis in tissue sections (MILAN technique) followed by image analysis, RT-PCR and shotgun proteomics. We observed that the brisk and non-brisk patterns are heterogeneous functional categories that can be further sub-classified into active, transitional or exhausted. The classification of primary melanomas based on the functional paradigm also shows correlation with spontaneous regression, and an improved prognostic value when compared to that of the brisk classification. Finally, the main inflammatory cell subpopulations that are present in the microenvironment associated with activation and exhaustion and their spatial relationships are described using neighbourhood analysis.

cisTopic: cis-regulatory topic modeling on single-cell ATAC-seq data.

We present cisTopic, a probabilistic framework used to simultaneously discover coaccessible enhancers and stable cell states from sparse single-cell epigenomics data ( http://github.com/aertslab/cistopic ). Using a compendium of single-cell ATAC-seq datasets from differentiating hematopoietic cells, brain and transcription factor perturbations, we demonstrate that topic modeling can be exploited for robust identification of cell types, enhancers and relevant transcription factors. cisTopic provides insight into the mechanisms underlying regulatory heterogeneity in cell populations.

Yin Yang 1 Orchestrates a Metabolic Program Required for Both Neural Crest Development and Melanoma Formation.

Increasing evidence suggests that cancer cells highjack developmental programs for disease initiation and progression. Melanoma arises from melanocytes that originate during development from neural crest stem cells (NCSCs). Here, we identified the transcription factor Yin Yang 1 (Yy1) as an NCSCs regulator. Conditional deletion of Yy1 in NCSCs resulted in stage-dependent hypoplasia of all major neural crest derivatives due to decreased proliferation and increased cell death. Moreover, conditional ablation of one Yy1 allele in a melanoma mouse model prevented tumorigenesis, indicating a particular susceptibility of melanoma cells to reduced Yy1 levels. Combined RNA sequencing (RNA-seq), chromatin immunoprecipitation (ChIP)-seq, and untargeted metabolomics demonstrated that YY1 governs multiple metabolic pathways and protein synthesis in both NCSCs and melanoma. In addition to directly regulating a metabolic gene set, YY1 can act upstream of MITF/c-MYC as part of a gene regulatory network controlling metabolism. Thus, both NCSC development and melanoma formation depend on an intricate YY1-controlled metabolic program.

The transcription factor Grainy head primes epithelial enhancers for spatiotemporal activation by displacing nucleosomes.

Transcriptional enhancers function as docking platforms for combinations of transcription factors (TFs) to control gene expression. How enhancer sequences determine nucleosome occupancy, TF recruitment and transcriptional activation in vivo remains unclear. Using ATAC-seq across a panel of Drosophila inbred strains, we found that SNPs affecting binding sites of the TF Grainy head (Grh) causally determine the accessibility of epithelial enhancers. We show that deletion and ectopic expression of Grh cause loss and gain of DNA accessibility, respectively. However, although Grh binding is necessary for enhancer accessibility, it is insufficient to activate enhancers. Finally, we show that human Grh homologs-GRHL1, GRHL2 and GRHL3-function similarly. We conclude that Grh binding is necessary and sufficient for the opening of epithelial enhancers but not for their activation. Our data support a model positing that complex spatiotemporal expression patterns are controlled by regulatory hierarchies in which pioneer factors, such as Grh, establish tissue-specific accessible chromatin landscapes upon which other factors can act.

SCENIC: single-cell regulatory network inference and clustering.

We present SCENIC, a computational method for simultaneous gene regulatory network reconstruction and cell-state identification from single-cell RNA-seq data (http://scenic.aertslab.org). On a compendium of single-cell data from tumors and brain, we demonstrate that cis-regulatory analysis can be exploited to guide the identification of transcription factors and cell states. SCENIC provides critical biological insights into the mechanisms driving cellular heterogeneity.

Comprehensive DNA methylation study identifies novel progression-related and prognostic markers for cutaneous melanoma.

BACKGROUND: Cutaneous melanoma is the deadliest skin cancer, with an increasing incidence and mortality rate. Currently, staging of patients with primary melanoma is performed using histological biomarkers such as tumor thickness and ulceration. As disruption of the epigenomic landscape is recognized as a widespread feature inherent in tumor development and progression, we aimed to identify novel biomarkers providing additional clinical information over current factors using unbiased genome-wide DNA methylation analyses. METHODS: We performed a comprehensive DNA methylation analysis during all progression stages of melanoma using Infinium HumanMethylation450 BeadChips on a discovery cohort of benign nevi (n = 14) and malignant melanoma from both primary (n = 33) and metastatic (n = 28) sites, integrating the DNA methylome with gene expression data. We validated the discovered biomarkers in three independent validation cohorts by pyrosequencing and immunohistochemistry. RESULTS: We identified and validated biomarkers for, and pathways involved in, melanoma development (e.g., HOXA9 DNA methylation) and tumor progression (e.g., TBC1D16 DNA methylation). In addition, we determined a prognostic signature with potential clinical applicability and validated PON3 DNA methylation and OVOL1 protein expression as biomarkers with prognostic information independent of tumor thickness and ulceration. CONCLUSIONS: Our data underscores the importance of epigenomic regulation in triggering metastatic dissemination through the inactivation of central cancer-related pathways. Inactivation of cell-adhesion and differentiation unleashes dissemination, and subsequent activation of inflammatory and immune system programs impairs anti-tumoral defense pathways. Moreover, we identify several markers of tumor development and progression previously unrelated to melanoma, and determined a prognostic signature with potential clinical utility.

Comparative oncogenomics identifies tyrosine kinase FES as a tumor suppressor in melanoma.

Identification and functional validation of oncogenic drivers are essential steps toward advancing cancer precision medicine. Here, we have presented a comprehensive analysis of the somatic genomic landscape of the widely used BRAFV600E- and NRASQ61K-driven mouse models of melanoma. By integrating the data with publically available genomic, epigenomic, and transcriptomic information from human clinical samples, we confirmed the importance of several genes and pathways previously implicated in human melanoma, including the tumor-suppressor genes phosphatase and tensin homolog (PTEN), cyclin dependent kinase inhibitor 2A (CDKN2A), LKB1, and others. Importantly, this approach also identified additional putative melanoma drivers with prognostic and therapeutic relevance. Surprisingly, one of these genes encodes the tyrosine kinase FES. Whereas FES is highly expressed in normal human melanocytes, FES expression is strongly decreased in over 30% of human melanomas. This downregulation correlates with poor overall survival. Correspondingly, engineered deletion of Fes accelerated tumor progression in a BRAFV600E-driven mouse model of melanoma. Together, these data implicate FES as a driver of melanoma progression and demonstrate the potential of cross-species oncogenomic approaches combined with mouse modeling to uncover impactful mutations and oncogenic driver alleles with clinical importance in the treatment of human cancer.

The effect of L-PRF membranes on bone healing in rabbit tibiae bone defects: micro-CT and biomarker results.

More insight into the biological fundamentals of leukocyte platelet-rich fibrin (L-PRF) guided healing is necessary to recommend its application, in particular in deficient bone sites that need to support implants. This study investigated the short-term bone healing effect of L-PRF treatment in cylindrical non-critical sized bone defects with 3 mm diameter and 6 mm depth in tibiae of 18 adult male New Zealand White rabbits. After a randomization process, 96 bone defects were prepared and half of them were filled with a L-PRF membrane, while untreated defects in the opposite tibia served as control group. The rabbits were euthanized after 7, 14 or 28 days of healing. The bone healing of the cortical and medullary areas was investigated by micro-CT, while the expression of molecular markers (RUNX2, VEGFA, COL1A2 and BMP2) was assessed by qRT-PCR. Treatment with L-PRF did not affect the micro-structural bone characteristics of the repaired bone tissue, except for a decrease in the trabecular connectivity at the cortical level after 14 days of healing. At this time, RUNX2 and VEGFA mRNA levels were significantly lower in the treated defects. L-PRF membranes thus had a temporary negative influence on the bone microarchitecture (Tb.Pf) and on the RUNX2 and VEGFA expression during early bone healing. Overall, L-PRF treatment did not enhance bone regeneration in these non-critical size defects after 28 days.