Multi-site tumor sampling highlights molecular
intra-tumor heterogeneity in malignant pleural
mesothelioma

Clément Meiller, François Montagne, Theo Hirsch, Stefano Caruso, Julien de Wolf, Quentin Bayard, Jean-Baptiste Assié, Léa Meunier, Yuna Blum, Lisa Quetel, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Clément Meiller, François Montagne, Theo Hirsch, Stefano Caruso, Julien de Wolf, et al. Multi-site tumor sampling highlights molecular intra-tumor heterogeneity in malignant pleural mesothelioma. Genome Medicine, 2021, 13 (1), pp.113. 10.1186/s13073-021-00931-w. inserm-03289952

HAL Id: inserm-03289952

https://inserm.hal.science/inserm-03289952

Submitted on 19 Jul 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# Multi-site tumor sampling highlights molecular intra-tumor heterogeneity in malignant pleural mesothelioma

Clément Meiller¹, François Montagne¹,², Theo Z. Hirsch¹, Stefano Caruso¹, Julien de Wolf¹,³, Quentin Bayard¹, Jean-Baptiste Assié¹,⁴,⁵, Léa Meunier¹, Yuna Blum⁶,⁷, Lisa Quetel¹, Laure Gibault⁸,⁹, Ecaterina Pintilie¹⁰, Cécile Badoual⁸,⁹, Sarah Humez¹¹,¹², Françoise Galateau-Sallé¹³, Marie-Christine Copin¹¹,¹⁴, Eric Letouzé¹, Arnaud Scherpereel¹⁵,¹⁶, Jessica Zucman-Rossi¹,⁸, Françoise Le Pimpec-Barthes¹,⁸,¹⁷, Marie-Claude Jaurand¹ and Didier Jean¹*

## Abstract

**Background:** Malignant pleural mesothelioma (MPM) is a heterogeneous cancer. Better knowledge of molecular and cellular intra-tumor heterogeneity throughout the thoracic cavity is required to develop efficient therapies. This study focuses on molecular intra-tumor heterogeneity using the largest series to date in MPM and is the first to report on the multi-omics profiling of a substantial series of multi-site tumor samples.

**Methods:** Intra-tumor heterogeneity was investigated in 16 patients from whom biopsies were taken at distinct anatomical sites. The paired biopsies collected from apex, side wall, costo-diaphragmatic, or highest metabolic sites as well as 5 derived cell lines were screened using targeted sequencing. Whole exome sequencing, RNA sequencing, and DNA methylation were performed on a subset of the cohort for deep characterization. Molecular classification, recently defined histo-molecular gradients, and cell populations of the tumor microenvironment were assessed.

**Results:** Sequencing analysis identified heterogeneous variants notably in NF2, a key tumor suppressor gene of mesothelial carcinogenesis. Subclonal tumor populations were shared among paired biopsies, suggesting a polyclonal dissemination of the tumor. Transcriptome analysis highlighted dysregulation of cell adhesion and extracellular matrix pathways, linked to changes in histo-molecular gradient proportions between anatomic sites. Methylome analysis revealed the contribution of epigenetic mechanisms in two patients. Finally, significant changes in the expression of immune mediators and genes related to immunological synapse, as well as differential infiltration of immune populations in the tumor environment, were observed and led to a switch from a hot to a cold immune profile in three patients.

<footnote>* Correspondence: didier.jean@inserm.fr
¹Centre de Recherche des Cordeliers, Inserm UMRS-1138, Sorbonne Université, Université de Paris, Functional Genomics of Solid Tumors, Paris, France
Full list of author information is available at the end of the article</footnote>
**Conclusions:** This comprehensive analysis reveals patient-dependent spatial intra-tumor heterogeneity at the genetic, transcriptomic, and epigenetic levels and in the immune landscape of the tumor microenvironment. Results support the need for multi-sampling for the implementation of molecular-based precision medicine.

**Keywords:** Thoracic tumor, Spatial molecular intra-tumor heterogeneity, Clonality, NF2 subclonal mutation, Tumor microenvironment

## Background

Malignant pleural mesothelioma (MPM) is a rare and highly aggressive tumor arising in the thoracic cavity. Exposure to asbestos is the main risk factor, and despite the ban of this mineral fiber in several countries, MPM remains a major public health problem worldwide. In most patients, MPM is an incurable cancer with a very poor prognosis, notably due to the ineffectiveness of conventional anti-tumor treatments. The reference treatment is based on systemic platinum-based chemotherapy combined with pemetrexed, a treatment that improves survival by only a few months even with the recent addition of the anti-VEGF therapy bevacizumab [1]. More recently, immunotherapy based on immune checkpoint inhibitors has shown survival benefits in particular in patients with non-epithelioid histology, without an accurate predictive biomarker for treatment response besides age and histology [2–4]. Despite this recent therapeutic progress, there is still an urgent need to develop a precision medicine approach taking into account MPM heterogeneity.

Like most solid tumors, MPM is a heterogeneous cancer with high variability between patients [5]. The heterogeneity of MPM was first described at the histologic level by defining three main types: epithelioid, sarcomatoid, and biphasic. Different histologic subtypes have been described especially for the epithelioid type including, but, not limited to, tubulopapillary, acinar, trabecular, solid, and micropapillary architectural subtypes [6, 7]. Interestingly, studies show that histologic subtype and grade of epithelioid MPM have a prognostic impact [6, 8]. More recently, the rare transitional type, which could also be a subtype of the sarcomatoid type, was described [9]. Large-scale omics and next-generation sequencing (NGS) studies have demonstrated MPM inter-tumor heterogeneity at the molecular level and led to molecular classifications into two to four subtypes [10–13]. More recent publications have shown that MPM heterogeneity is well-described by a continuum [10, 14]. Recently, we defined histo-molecular gradients, which also takes into consideration intra-tumor heterogeneity, by determining the proportions of epithelioid-like and sarcomatoid-like cellular entities (E.score and S.score), related to the two extreme histological types of MPM, within tumor samples. These gradients have a high prognostic value and are of interest to guide

therapeutic strategies, including targeted therapies and immunotherapies [10].

MPM is a solid tumor characterized by a diffuse locoregional growth within the pleural cavity. A polyclonal origin has been described and subclonal cell populations with specific mutations have been evidenced [15, 16]. So far, only two studies have focused on spatial genetic heterogeneity [17, 18]. Variability in the mutational load without involving key genes of mesothelial carcinogenesis was first observed between anatomical locations in a series of six patients [17]. Very recently, a series of nine patients also highlighted heterogeneous mutations between anatomical locations [18]. Some evidence supports spatial heterogeneity concerning the tumor microenvironment [5]. Both the above studies showed distinct T-cell repertoires depending on the tumor anatomic region [17, 18]. However, contrary to inter-tumor microenvironmental heterogeneity in terms of stromal and immune cell infiltration [10, 19], intra-tumor cell heterogeneity of the tumor microenvironment has not yet been reported.

In the present study, we have gone further in the characterization of MPM molecular spatial heterogeneity in the largest series of patients (16 cases) studied so far. We not only defined the heterogeneity at the genetic level and, for the first time, characterized the spatial dysregulation of gene expression, epigenetic changes, and tumor microenvironment differential infiltration.

## Methods

### Patients

Tumor samples were collected from up to four distinct anatomical sites (apex, side wall, costo-diaphragmatic, and highest metabolic sites detected by positron emission tomography (PET) scan when present, shown in Additional file 1: Figure S1) in a series of 16 patients who had diagnostic biopsies or surgery resections for MPM in two French hospitals (CHRU of Lille and Hôpital Européen Georges Pompidou of Paris). For each patient, the tumor location is indicated (A: apex; B: side wall; C: costo-diaphragmatic; D: highest metabolic site). Non-tumoral samples (blood or muscle of the chest wall) were available for 9 patients. Patients were diagnosed between 2014 and 2017 and tumors were certified by the French National Pathology Expertise Network (Mesopath) as MPM [20]. Histologic type, subtype, and
grade were determined by MesoPath expert pathologists and part of the series was reviewed according to the WHO 2021 update [21]. The experiments were undertaken with the understanding and written consent of each subject. The study methodologies were conformed to the standards set by the Declaration of Helsinki and approved by a local medical ethics committee (CPP Ile-de-France II). The sampling procedures were approved by the French research ministry (CODECOH no. DC-2016-2771). Samples were annotated with detailed clinico-pathological and epidemiologic information obtained from pathology reports (Additional file 1: Table S1). Based on tumor purity and a quality check of extracted nucleic acids, two tumor samples per patient were used for further analysis. Three metrics were used to evaluate tumor purity: (i) pathologist assessment of tumor cellularity, (ii) prediction by RT-qPCR of non-tumor components using WISP deconvolution method [10], and (iii) frequencies of variants determined by targeted NGS. For the few patients with more than 2 qualified tumor samples, paired samples with the closest tumor purity were kept. In addition, MPM primary cell lines were generated from five patients in our laboratory, from fresh tumor samples collected at sites distinct from the four previous locations. They were established, cultured as previously described [22], and used at low-passage numbers (6 to 10 passages). All tumor samples and cell lines were analyzed by targeted NGS, 9 paired tumor samples by whole exome sequencing (T199LE-A/B, T200LE-A/C, T201LE-A/B, T203LE-A/C, T225LE-A/D, T227LE-A/D, T277HP-A/C, T278HP-A/C, T333HP-A/C), 8 by RNA sequencing (T199LE-A/B, T201LE-A/B, T203LE-A/C, T225LE-A/D, T227LE-A/D, T277HP-A/C, T278HP-A/C, T333HP-A/C), and 5 by methylation profiling (T203LE-A/C, T227LE-A/D, T277HP-A/C, T278HP-A/C, T333HP-A/C).

## Nucleic acid extraction

For tumor samples, a preliminary step of tissue disruption and cell lysis was achieved using TissueLyser II (Qiagen, Courtaboeuf, France). Genomic DNA and total RNA from tumor samples and cell pellets, obtained from cultures of MPM primary cell lines at passages 6 to 10, were extracted using the AllPrep DNA/RNA/miRNA Universal kit (Qiagen) according to the manufacturer's protocol. For the fresh tumor sample used to establish MPM_83, automated DNA extraction was performed on multiple sections (n = 34) of the frozen preserved sample following the protocol of the Maxwell 16 Tissue DNA purification kit and the Maxwell instrument (Promega, Charbonnières-les-Bains, France). DNA and RNA quantifications were done by fluorescence measurements (Hoechst dye) on a FLUOstar Omega microplate reader (BMG Labtech, Champigny sur Marne, France) and by

spectrometry on a NanoDrop-1000 (Ozyme, Saint-Cyr l'Ecole, France), respectively. We used agarose gel migration for DNA to ensure the absence of excessive degradation before sequencing and methylation analysis. We assessed RNA integrity on a Fragment Analyzer (Agilent Technologies, Courtaboeuf, France) before RNA-seq analysis.

## Targeted next-generation sequencing (NGS) and variant calling

We performed targeted NGS using our in-house protocol recently published in detail to sequence 21 genes and the TERT promoter on the whole series of samples (16 patients) [22, 23]. Briefly, library preparation was based on a multiplex PCR enrichment and sequencing was achieved by a MiSeq instrument (Illumina, Evry, France). FASTQ files were generated by the Illumina MiSeq Reporter software. Primer sequences were removed using the fastx_trimmer function from the fastx toolkit. Reads were aligned on the human genome assembly hg19/GRCh37 using the Burrows-Wheeler aligner (BWA) and bam files were generated using samtools. Variant calling was performed with Unified Genotyper and variant annotation was obtained using the Oncotator annotation algorithm, along with the ensembl Variant Effect Predictor (VEP) algorithm and the Annovar annotation. Finally, based on these annotations, all somatic variants with a functional consequence were checked by visualization using Integrative Genomic Viewer (IGV) software (Broad Institute, Cambridge, MA, USA). Variants of interest were validated either using whole exome sequencing (WES) data or by a second independent targeted NGS. Genome coordinates were converted to human genome assembly hg38/GRCh38 using the UCSC LiftOver online tool [24].

## Whole exome sequencing (WES) and variant calling

Library preparation, exome capture, sequencing, and data analysis were performed by IntegraGen (Evry, France). We screened all paired biopsies with non-tumoral samples available (9 patients) [23]. Briefly, genomic DNA was captured either using the Agilent in-solution enrichment methodology with their biotinylated oligonucleotide probe library (SureSelect Clinical Research Exome V2, Agilent Technologies) or the Twist Human Core Exome Enrichment System (Twist Bioscience, San Francisco, USA), followed by paired-end 75 base massively parallel sequencing on Illumina HiSeq4000. Sequence capture, enrichment, and elution were performed according to the manufacturer's instructions and protocols without modification, except for the library prepared using NEBNext® Ultra II kit (New England Biolabs, Evry, France). Image analysis and base calling were performed using Illumina Real-Time Analysis
(2.7.7) with default parameters. Sequence reads were mapped to the human genome build hg38/GRCh38 by using the BWA tool. The duplicated reads were removed using sambamba tools. Variant calling, allowing the identification of genetic alterations, as well as SNV (single nucleotide variation) and small insertions/deletions (up to 20 bp), was performed via the Broad Institute's GATK Haplotype Caller GVCF tool (3.7) for constitutional DNA and via the Broad Institute's MuTect tool for somatic DNA. Ensembl's VEP (variant effect predictor) program processed variants for further annotation [25, 26]. An in-house post-processing workflow was applied to filter out candidate somatic mutations that were more consistent with artifacts or germline mutations. Finally, the variants were validated by visualization using IGV software. Two bioinformatics predictions for missense pathogenicity were used: SIFT (5.2.2) and PolyPhen (2.2.2) and damaging variants were considered if it was predicted by at least one tool. The circular binary segmentation algorithm implemented in the Bioconductor package DNAcopy (DNAcopy 1.32.0) as well as FACET R package (v.0.6.1.) were used to reconstruct copy-number profiles from WES data [27, 28]. For FACET analysis, single nucleotide polymorphisms (SNP) count matrix for both tumoral and non-tumoral samples was obtained by processing bam files with snp-pileup (arguments: -q15; -Q20; -P100; -r20,0). Then, SNP matrix was processed using preProcSample function (cval = 25, snp.nbhd=250, ndepth=30) to generate log-R-ratio and the segmentation. ProcSample (cval=150, min.-nhet=5) was then used to estimate the wild-type 2-copy state. Finally, emcncf (min.nhet=5) function was used to estimate sample ploidy and purity. R package ggplot was used to make pangenomic graphical representations. Cancer-related genes were based on Tier 1 of the Cancer Gene Census (CGC) database [29, 30].

## Clonality prediction

For each mutation, we used the Palimpsest R package [31] to estimate the fraction of tumor cells harboring this variant (cancer cell fraction, CCF) in each tumor, as previously described [32]. In addition, the emcncf (min.-nhet=5) function of FACET package (v.0.6.1.) was used to determine cellular fraction and allele-specific copy-numbers per segments from CNV data [28]. T200LE was excluded from this analysis as results were not reliable probably due to low tumor purity. CCF for CNV was calculated by dividing a cellular fraction by tumor purity. For pairs of tumor samples from a given patient, the comparison of the CCF distribution between samples was used to sort mutations or CNV among clonal shared, subclonal shared, clonal private, and subclonal private variants. Variants or CNV were considered subclonal if the CCF value of both tumor samples were

lower than 0.5 and, for variants only, if the upper limit of the confidence interval was below 1. Shared variants were found in both paired tumor samples in contrast to private variants.

## RNA sequencing

Library preparation and sequencing were performed by IntegraGen on paired biopsies of 8 patients [23]. Libraries were prepared using the TruSeq Stranded mRNA kit (Illumina) or the NEBNext Ultra II Directional RNA Library Prep Kit (New England Biolabs), according to the supplier's recommendations. Briefly, the key steps of both protocols were (i) purification of PolyA-containing mRNA molecules using poly-T oligo attached magnetic beads from 1 µg total RNA, (ii) RNA fragmentation using divalent cations at an elevated temperature to obtain approximately 300 bp fragments, (iii) double-strand cDNA synthesis, and (iv) Illumina adapter ligation and cDNA library amplification by PCR for sequencing. Paired-end 75 base massively parallel sequencing was then carried out on an Illumina HiSeq4000. FASTQ files were aligned to the human reference genome hg38/GRCh38 using TopHat2 V.2.0.14 [33]. We removed reads mapping to multiple locations, and we used HTSeq [34] to obtain the number of reads associated with each gene in the Gencode v27 database. The Bioconductor DESeq2 package [35] was used to import raw HTSeq counts for each sample into R statistical software (R Foundation for Statistical Computing, Vienna, Austria) and to apply variance stabilizing transformation to the raw count matrix. FPKM values (number of fragments per kilobase of exon model and millions of mapped reads) were calculated by normalizing the count matrix for the library size estimated with the DESeq2 package and the coding length of each gene. Differentially expressed genes were determined using both FPKM and DESeq2-normalized RNA-seq data: genes that displayed a FPKM value below 1 in both paired biopsies were excluded and we considered as differentially expressed only genes with a difference of DESeq2 expression values between paired samples above 1. Hierarchical clustering was done using cosine distance and Ward's linkage method in R statistical software with the 500 most differentially expressed genes. Network analysis of the distribution of the differentially expressed genes among tumor samples was performed using the web-based interactive tool DiVenn [36].

## Gene fusion detection

Fusions detected by TopHat2 (--fusion-search --fusion-min-dist2000 --fusion-anchor-length13 --fusion-ignore-chromosomes chrM) were filtered using the TopHatFusion pipeline [37]. We discarded artifacts based on fusion redundancy among samples from
different patients. Fusions validated by BLAST and
with at least 20 split-reads or pairs of reads spanning
the fusion event in at least one of the paired tumor
samples were retained. Shashimi plots were generated
using junctions-bed provided by TopHat2 to deter-
mine the junctions supported by splitted reads and
bedtools genomcov function (with -bga -split op-
tions) to calculate the coverage from the aligned bam
files.

Molecular classification and histo-molecular score
predictions

Expression data from RNA sequencing (DESeq2-normalized counts) were used to predict molecular subgroups and histo-molecular scores (E/S.score) based on cen- troids and deconvolution approaches, respectively, as de- scribed previously [10].

Methylation profiling

Experimental works and data analysis were performed by IntegraGen on paired biopsies of 5 patients [38]. First, bisulfite conversion of 500 ng of DNA samples was performed with the Zymo EZ DNA Methylation Kit (Ozyme) according to the manufacturer's instructions. Then, methylation analysis was performed on the Infinium MethylationEPIC BeadChip Kit (Illumina) following the supplier's recommendations. Briefly, samples were denatured, neutralized, and incubated in the Illumina hybridization oven for 20-24 h at 37°C to uniformly amplify genomic DNA, then enzymatically fragmented. Following the fragmentation step, samples were hybridized to the BeadChip and incubated at 48°C for 16-24 h in the Illumina hybridization oven. Finally, non-hybridized and unspecific hybridized DNA samples were washed from the BeadChip, and labeled nucleotides were added to extend and stain primers hybridized to the samples. The Illumina iScan system (with the iScan Control Software) was used to scan Illumina MethEpic BeadChips and raw data were analyzed with GenomeStudio software v2011 according to Illumina's recommendations to generate result files with beta-values. Loci were considered as differentially methylated between paired tumor samples when both detection *p*-values were lower than 0.05 and a difference of at least 0.2 was observed between beta-values of each tumor sample. To compare with genes previously reported to be correlated to the E.score or the S.score [10], we only considered genes with expression changes in the expected sense according to the histo-molecular scores of the paired biopsies. Hierarchical clustering was done using cosine distance and Ward's linkage method in R statistical software with the 500 CpG showing the most variability in their methylation beta-value.

Pathway dysregulation analysis

Signal pathway dysregulation analysis was performed by two complementary approaches : over-representation analysis and single-sample Gene Set Enrichment Analysis (ssGSEA). Differentially expressed or methylated genes lists obtained for each patient were used as input for over-representation analysis in the web-based Gene Set Analysis Toolkit using the KEGG, Reactome, or GeneOntology Biological Process non-redundant databases [39]. We only considered pathways with false discovery rate (FDR) adjusted *p*-values lower than 0.05. Gene ratio is the proportion of differentially expressed genes regarding the total number of genes included in a given pathway. SsGSEA were calculated as gene-set variation analysis (GSVA) enrichment scores using the DESeq2-normalized RNA-seq data, after excluding genes with a FPKM value below 1 for all samples, and the GSVA package [40]. The absolute value of the delta of each ssGSEA score between paired tumor samples (ssGSEA_score) was used to assess intra-tumor heterogeneity.

Tumor microenvironment estimation

The microenvironment cell population counter (MCP-
counter) method [41] was used to compute scores of in-
filtration for different cell populations from DESeq2-
normalized RNA-seq data. To optimize the prediction of
the relative abundance of stromal and immune cells
within MPM tumors, the genes included in the signature
of each cell population were selected based on their ab-
sence of expression in 22 MPM cell lines. IGKC and
CHRM3-AS2 were also excluded as gene expression was
not present in all datasets. The complete list of genes
used for each population is available in Additional file 1:
Table S2. For integrative analysis with public datasets, a
total of 295 samples from three different series of RNA-
seq were combined in the integrated analysis, including
the 16 paired tumor samples, 209 samples from Bueno
et al. [11, 42], and 70 samples from TCGA [13] down-
loaded through the Broad Institute TCGA GDAC fire-
hose tool [43]. Then, we standardized gene expression
separately to have a mean of 0 and a standard deviation
of 1 per gene in each dataset. Statistical analysis and data
visualization were performed using R software. Unsuper-
vised hierarchical clustering was performed using cosine
distance and Ward's linkage method. The Wilcoxon test
was used to estimate the difference in the abundance of
immune populations within paired tumor samples.

Immunohistochemical staining

Formalin-fixed, paraffin-embedded (FFPE) tumor biopsies (3 biopsies at distant sites per patient) were sectioned at a thickness of 3 µm and stained on positively charged glass slides. Deparaffinization, rehydration, and
antigen retrieval were performed by CC1 (prediluted; pH 8.0) antigen retrieval solution (Ventana Medical Systems, Inc.) on the Ventana BenchMark ULTRA automated slide stainer for 32 min at 100 °C. Specimens were incubated with primary antibodies anti-CD3 (clone LN10; Leica; dilution 1:100), anti-CD8 (clone C8/144B; Dako, dilution 1:25), and anti-CD20 (clone L26; Diagnostics; dilution 1:150) followed by visualization with the Ultraview DAB IHC Detection Kit. The specimens were then counterstained with hematoxylin II and bluing reagent (Ventana) and coverslipped. Each IHC run contained a positive control.

# Results

## Mutational intra-tumor heterogeneity

Paired tumor samples collected prior to chemotherapy treatment at two distinct anatomical sites from 16 MPM patients and five MPM-derived primary cell lines were analyzed by different NGS and omics methods (Fig. 1a). Patient characteristics are detailed in Additional file 1: Table S1. In summary, patients were chemo naïve, mostly male (81%), had a past exposure to asbestos (62%), and were diagnosed at a median age of 74 years. One tumor was biphasic and the others were epithelioid, the most frequent histologic type. Even if the epithelioid tumors belong to different subtypes and grades, our series do not fully reflect the histologic diversity of MPM. Germ-line mutations in cancer-related genes for 9 patients are reported in Additional file 1: Table S3A. Among genes previously reported as altered by pathogenic or likely pathogenic germline mutations in MPM [44], damaging variants were found in *SDHA* and *PALB2* genes in patients T225LE and T277HP, respectively.

To investigate mutational intra-tumor heterogeneity, we performed targeted NGS on a panel of key genes of mesothelial carcinogenesis [22] for the entire series and WES for patients with non-tumoral DNA available to go deeper in their genetic characterization. Most of the protein-altering somatic variants, validated using multiple sequencing approaches and IGV visualization, were common between paired samples from a given patient (Fig. 1b, c; Additional file 1: Table S3B-D). However, we found intra-tumor heterogeneous variants in the well-known mesothelioma driver gene *NF2*. Among eight patients harboring *NF2* mutations, three displayed intra-tumor heterogeneity. The first patient (T227LE) had a *NF2* nonsense variant in the derived primary cell line (MPM_83), although it was not called in the paired biopsies. Therefore, we investigated the presence of the variant in the fresh tumor sample used to establish MPM_83 by sequencing multiple sections of this sample and found the *NF2* variant at different frequencies, confirming the subclonal status of the *NF2* mutation (Additional file 2: Figure S2, Additional file 1: Table S3E). For

the second patient (T278HP), we detected an *NF2* splice site variant in the sample localized at the apex but it was absent in the paired costo-diaphragmatic sample and in the derived primary cell line (MPM_70). Finally, the third patient (T333HP) presented an *NF2* nonsense variant at low frequency in both biopsies, but the mutation was not found in the derived primary cell line (MPM_71), suggesting that this mutation was subclonal and probably absent in the fresh tumor sample used to establish MPM_71. In addition to *NF2*, we also detected intra-tumor heterogeneous mutations in three cancer-related genes [29]: (i) a frameshift insertion variant in *CTNNB1* for patient T225LE, (ii) a missense variant in the protein tyrosine phosphatase receptor *PTPRT* for patient T227LE (however, based on RNA-seq data, this gene does not seem to be expressed in both samples of this tumor), and (iii) a nonsense variant in the transcriptional coactivator *PSIP1* for patient T333HP, to our knowledge the first mutation in this gene described for MPM (Fig. 1b).

The clonality of the protein-altering variants was also assessed to investigate subclonal tumor populations that could be in different proportions among paired samples. To perform a robust analysis of clonality, we calculated the cancer cell fraction of each variant, by normalizing variant frequencies with copy number and tumor purity (Additional file 1: Table S4A). As expected, we found most of the cancer-related gene variants to be clonal and shared by both paired samples except for the four genes mentioned previously (Additional file 2: Figure S3). Tumor clonality was heterogeneous between patients, who showed a variable number of private variants ranging from 0 to 12 (Additional file 1: Table S4B). Strikingly, we found a subclonal tumor population present in both biopsies in patient T333HP supported by six variants, including *NF2* (Fig. 1d). Of note, the variant frequencies in the clonal tumor population (*ARID2*, *SETD2*, and *TERT* promoter) and in the subclonal population (*NF2*) were validated by targeted NGS for this patient (Additional file 1: Table S3B). This suggests the polyclonal dissemination of these tumors throughout the thoracic cavity. Furthermore, shared subclonal variants were also identified in three other patients, i.e., T199LE, T203LE, and T225LE (Additional file 2: Figure S3).

## Copy number variations and fusion transcripts

Because MPM is often characterized by chromosomal alterations, we took advantage of WES coverage and RNA-seq data to analyze copy number variations (CNV) and fusion transcripts, respectively. First, we found that pangenomic CNV between paired samples were very similar using both DNA copy or FACET tools, notably in the *CDKN2A*, *BAP1*, and *NF2* loci (Additional file 2: Figures S4 and S5). Segmentation and CNV clonality
<figure><img src="image_4.png" /></figure>

provided by FACET analysis identified mainly clonal
shared CNV, but also some clonal private (n = 22), sub-
clonal private (n = 4), and subclonal shared (n = 1) CNV

(Additional file 1: Table S4C). However, all homozygous
deletions were clonal and shared by paired tumor sam-
ples, and, as expected, often located in the CDKN2A
locus. Interestingly, we identified a CNV with a subclonal status in both paired tumor biopsies in patient T333HP, supporting previous results based on variant clonality analysis. Then, we searched for fusion transcripts in RNA-seq data that we generated for the same series of samples screened with WES, except for tumor samples from patient T200LE due to poor RNA quality in addition to very low tumor content. In three patients, we detected fusions involving a cancer-related gene, such as *STK11* and *KDM6A*, quoted as altered in MPM in the Cosmic database (release v91), or *TFG* reported in fusions in other cancers (Additional file 1: Table S5, Additional file 2: Figure S6A) [45]. Interestingly, we found that the fusion transcript in patient T277HP involving *TFG* and *ADGRG7* led to abnormal gene expression of the adhesion G-protein coupled receptor *ADGRG7*, a gene not commonly expressed in MPM (Additional file 2: Figure S6B). For these three fusions, both paired samples harbored the rearrangement. To note, we also detected four fusion transcripts restricted to one of the paired samples in three patients, but they did not involve cancer-related genes. In conclusion, we did not observe structural intra-tumor heterogeneity at the genomic level with an evident contribution to cancer development.

## Differential gene expression and signal pathway dysregulation

Gene expression analysis, from the RNA-seq data of paired samples from eight patients, indicated global gene expression as being patient-specific, as shown by unsupervised clustering (Additional file 2: Figure S7). The number of differentially expressed (DE) genes between paired samples ranged from 647 to 2119 genes, regardless of tumor cellularity (Spearman correlation test, R squared=0.01, *p* = 0.79), with a majority of protein coding genes (Additional file 2: Figure S8A). A high proportion of these protein coding genes was found to be differentially expressed in several tumor samples (Additional file 1: Table S6A, Additional file 2: Figure S8B). However, shared DE genes were distributed throughout all tumor samples (Additional file 2: Figure S8C). We performed signal pathway overrepresentation analysis using the KEGG, Reactome, and GeneOntology databases on protein coding genes differentially expressed between paired samples. Among the top recurrent pathways detected in all these databases, we found two main families, the first one related to cell adhesion and extracellular matrix organization and the second to immune communication (Additional file 1: Table S6B-D).

We also analyzed signal pathway dysregulation by ssGSEA across the whole cohort (Additional file 1: Table S6E). Recurrent major changes between paired tumor

samples concern immune pathways, more precisely related to immunological synapses involving specific immune cells infiltration such as antigen-presenting cells (dendritic and B cells) and T cells.

## Histo-molecular heterogeneity

We previously proposed classification in four molecular subtypes and histo-molecular gradients determining the proportion of epithelioid-like and sarcomatoid-like components (E/S.scores) to describe MPM inter- and intra-tumor heterogeneity [10]. Based on RNA-seq data, we first predicted the molecular subtypes and obtained similar subgroups between paired biopsies (Fig. 2a). Then, we estimated the E/S.scores and observed changes greater than 10% in the E/S.scores among paired samples in the three patients T199LE, T227LE, and T278HP. Interestingly, these three patients showed differential enrichment for the pathways associated with cell adhesion and extracellular matrix (Fig. 2a), which makes sense in the context of our previous publication [10]. As expected, patients T201LE, T203LE, T277HP, and T333HP displayed similar E/S.scores between paired samples and did not show major changes in these pathways. In contrast, patient T225LE showed a slight increase in the S.score, but a differential expression of genes enriched in these pathways. ssGSEA data also show a globally higher variation of these pathways between paired tumor samples of patients T225LE, T227LE, and T278HP and to a lesser extent T199LE compared to the four others (Additional file 2: Figure S9).

## Epigenetic intra-tumor heterogeneity

To evaluate intra-tumor heterogeneity at the epigenetic level, we performed methylome analysis in five patients with sufficient remaining quantities of DNA. As for the transcriptome analysis, the global methylation profile, illustrated by unsupervised clustering, remained very specific for each patient (Additional file 2: Figure S10A). Differential methylation analysis between paired samples showed that the majority of differentially methylated (DM) CpG were associated with protein coding gene loci in all patients (Additional file 1: Table S7A, Additional file 2: Figure S10B). However, this analysis highlighted that the number of DM CpG between paired samples varied by patient, with two patients (T227LE and T278HP) showing a high number of DM CpG, 5476 and 19178 respectively (Additional file 2: Figure S10C). Taking into account protein coding gene expression, a potential impact of methylation was only observed in these two patients, as DM CpG were identified in 18% and 41%, respectively, of DE genes (Fig. 2b). Interestingly, among the
A

<table><thead><tr><th>Patient</th><th>T199LE</th><th>T201LE</th><th>T203LE</th><th>T225LE</th><th>T227LE</th><th>T277HP</th><th>T278HP</th><th>T333HP</th></tr><tr><th>Localisation</th><th>A</th><th>B</th><th>A</th><th>B</th><th>A</th><th>C</th><th>A</th><th>C</th></tr></thead><tbody><tr><th>Classification (2 subtypes)</th><td>C1</td><td>C1</td><td>C1</td><td>C1</td><td>C1</td><td>C1</td><td>C1</td><td>C1</td></tr><tr><th>Classification (4 subtypes)</th><td>C1B</td><td>C1B</td><td>C1A</td><td>C1A</td><td>C1B</td><td>C1B</td><td>C1A</td><td>C1A</td></tr><tr><th>E.score</th><td>92%</td><td>78%</td><td>83%</td><td>82%</td><td>81%</td><td>83%</td><td>95%</td><td>91%</td></tr><tr><th>S.score</th><td>8%</td><td>22%</td><td>17%</td><td>18%</td><td>19%</td><td>17%</td><td>5%</td><td>9%</td></tr></tbody></table><table><thead><tr><th>Cell adhesion</th><th>Positive regulation of cell adhesion [GO:0045785]</th><th>Cell adhesion mediated by integrin [GO:0033627]</th><th>Integrin activation [GO:0033622]</th><th>Cell adhesion molecules (CAMs) [hsa04514]</th><th>Focal adhesion [hsa04510]</th><th>Integrin cell surface interactions [R-HSA-216083]</th></tr></thead><tbody><tr><th>Extracellular matrix</th><th>Extracellular structure organization [GO:0043062]</th><th>ECM-receptor interaction [hsa04512]</th><th>Extracellular matrix organization [R-HSA-1474244]</th><td></td><td></td><td></td></tr></tbody></table>

B

<figure><img src="image_5.png" /></figure>

C

<figure><img src="image_6.png" /><figcaption>**Fig. 2** Intra-tumor heterogeneity at the transcriptomic and epigenetic levels. **a** Molecular heterogeneity (molecular classifications and histo-molecular gradients) of the paired tumor samples and the over-represented pathways linked to cell adhesion and the extracellular matrix are shown in the table. For each patient, the tumor location is indicated (A: apex; B: side wall; C: costo-diaphragmatic; D: highest metabolic site). The molecular classifications into two to four subtypes and the E/S.scores were predicted based on RNA-seq data and are colored in blue or green depending on the subtypes and with a red gradient for the E/S.scores. Pathway over-representation is indicated by a circle with the size and a color proportional to the gene ratio and the FDR *p*-values, respectively. **b, c** Based on transcriptome and methylome analysis, the differentially expressed protein coding genes with an associated differentially methylated CpG (DE_DM genes) between paired tumor samples were determined. The percentage of DE_DM genes among all the differentially expressed protein coding genes is shown in the histogram for each patient. The number of DE_DM genes is indicated at the right of the histogram bars (**b**). The proportion of protein coding genes previously shown to be correlated to the E/S.scores [10] in all DE_DM genes is indicated in the pie charts for patients T227LE and T278HP (**c**). DE: differentially expressed; DM: differentially methylated</figcaption></figure>

DE protein coding genes associated to DM CpG, 18% and 47%, respectively, in patients T227LE and T278HP, are genes previously identified to be correlated with the E/S.scores of the histo-molecular gradients (Fig. 2c). Furthermore, in T278HP, these genes were enriched in previously identified pathways involved in cell adhesion and extracellular matrix organization (Additional file 1: Table S7B; Additional file 2: Figure S10D). These results support the notion that epigenetic regulation could contribute to the

variations occurring in the proportions of epithelioid-like and sarcomatoid-like components in patient T278HP and to a lesser extent in patient T227LE, in agreement with our previous studies demonstrating that histo-molecular gradients are related to epigenetic regulation [10].

## Tumor microenvironment

We next analyzed the immune landscape changes in paired samples. Comparison based on over-
representation analysis of the immune system-related pathways redundant among the KEGG, Reactome, and GeneOntology databases (Additional file 1: Table S6) highlights pathways belonging to immune cell communication such as cytokines and chemokines (Fig. 3a). The ssGSEA analysis points out changes in immunological synapse pathways (Fig. 3b). Overall, dysregulation in immune pathways between paired tumor samples was higher in T199LE, T225LE, T227LE, and to a lesser extent T201LE. Several immune checkpoints also showed differential expression between paired samples (Fig. 3c). In particular, a fold change greater than 2 was found in immune checkpoints targeted by immunotherapy, i.e., PDL1 (CD274) for T277HP and T278HP, CTLA4 (CTLA4) for T227LE, and PD1 (PDCD1) for T278HP. In addition, we estimated the relative proportions of immune and stromal populations using MCP-counter defined biomarkers (Additional file 1: Table S2) from RNA-seq data to determine the consistency of the tumor microenvironment between paired samples. Clustering analysis using the estimated cell populations from paired samples integrated with tumor samples from the Bueno and TCGA series [11, 13] separated tumor samples in two clusters, one with a high level of immune cell infiltration corresponding to a hot tumor immune profile, and the other with a lower level of infiltration corresponding to a cold profile (Fig. 4a). Three paired tumor samples from patients T199LE, T225LE, and T227LE, which displayed the highest dysregulation in immune pathways using both pathway analyses, were distributed separately between hot and cold tumors, indicating major changes in terms of the tumor microenvironment composition. Clustering of paired samples without the integration of other tumor samples also separated the tumor paired samples of patients T199LE, T225LE, and T227LE between the cold and hot phenotypes (Additional file 2: Figure S11). These paired samples showed significantly different relative proportions of stromal and immune cell populations with an increase in the infiltration of all cell populations in the hot immune profile paired sample compared to the cold immune profile sample (Fig. 4b). More complex changes with an increase or decrease in certain cell populations between paired samples were also observed in other patients (Additional file 2: Figure S12A). Individual comparison of cell populations showed that natural killer cells, cytotoxic lymphocytes, myeloid dendritic cells, T cells, B cells, and T CD8 cells were variable in paired tumors switching from a cold to a hot immune profile (Additional file 2: Figure S12B). SsGSEA analysis of the aggregated pathways linked to specific immune cell populations confirmed the greatest changes in the infiltration within paired tumors with a hot/cold mixed profile. Comparison between paired tumors with a hot/cold

mixed profile and the others displayed significant differences (Mann-Whitney test, *p* = 0.04) for natural killer cells, myeloid dendritic cells, and T cells (Additional file 2: Figure S13). IHC using anti-CD3, anti-CD8, and anti-CD20 antibodies on 3 FFPE tumor samples per patient, biopsied at distant anatomic sites, showed heterogeneous infiltration of T and B cells in T199LE patient and to a lesser extent in T225LE and T227LE patients (Additional file 2: Figure S14).

## Discussion

The present study provides a comprehensive description of spatial intra-tumor heterogeneity in MPM, along with new insights into tumor evolution and clues to the impact on therapy.

By a robust analysis of clonality restricted to validated somatic protein-altering variants, we confirmed the occurrence of MPM populations with different mutation profiles, depending on the anatomic site in the same patient, as previously observed [17, 18]. All tumors harbored many shared clonal mutations, supporting the notion that MPM was derived from a single cell in our series, contrary to the polyclonal origin previously described [15]. Along with these shared clonal variants, the presence of private clonal or subclonal variants indicates that tumor cells spread in the thoracic cavity and continue to evolve separately at different locations. One important finding was the existence of subclonal tumor populations shared between anatomic sites for at least one patient and likely in three others, which supports polyclonal dissemination throughout the thoracic cavity in some MPM patients. Evidence of polyclonal dissemination has been observed in metastases for several cancers [46]. This is possibly related to multiple waves of migrating cells or to the implantation of circulating tumor cell clusters composed of genetically distinct clones that have been described in biological fluids [47]. For MPM, numerous large tissue fragments and multicellular balls or berry-like clusters of cells are observed by cytology in the pleural fluid of patients [48, 49]. Our results suggest that these multicellular structures of tumor cells participate in the spread of MPM in the thoracic cavity.

NF2 is one of the main mutated tumor suppressor genes in MPM. As it plays a key role during mesothelial carcinogenesis, inactivating variants were expected to be clonal [50]. However, we clearly demonstrate the presence of subclonal NF2 mutations in three patients, characterized by clonal mutation in other key MPM mutated cancer genes including *BAP1*, *SETD2*, and the *TERT* promoter. The existence of MPM subclones harboring NF2 mutations was first suggested by a study identifying one cell line with homozygous NF2 mutation, although this was almost undetectable in the original tumor in a series of nine MPM patients [16]. Heterogenous NF2
<figure><img src="image_7.png" /></figure>

mutations between different tumor regions were also reported in two out of nine patients in another study, but the subclonality was not confirmed by a rigorous CCF analysis as we performed here [18]. Altogether, these

results suggest that subclonal *NF2* mutations are frequent in MPM and support that *NF2* inactivation could be a late event during MPM development. Merlin encoded by the *NF2* gene, is a multifunctional protein
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Differential infiltration of stromal and immune cell populations. a Unsupervised clustering of the paired tumor samples of eight patients with 209 and 70 tumor samples from the Bueno and TCGA series, respectively, was performed based on cell populations determined by the MCP-counter method. The paired tumor samples of the eight patients (series U1138) are indicated by a color code at the top of the heatmap as well as the series of each tumor sample. b The violin plots show the normalized MCP-counter values of immune and stromal cell populations between paired tumor samples for patients T199LE, T225LE, and T227LE, which are characterized by a hot/cold mixed immune profile. Each cell population is indicated in the box plots by a color point connected by a dotted line between paired tumor samples. The *p*-values of the Wilcoxon test comparing distribution between paired tumor samples are indicated at the top of the violin plot</figcaption></figure>

exhibiting well-established tumor-suppressive function through several cellular processes, not only linked to the activation of Hippo pathway, but also to the inhibition of PI3K/AKT/mTOR pathway as well as to regulatory

functions in the nucleus [50]. Although *NF2* inactivation is clearly a driver of mesothelial carcinogenesis, it is possible that this inactivation does not play a role in tumor initiation for MPM, but will promote the development
of a more aggressive tumor. This hypothesis is in line
with several findings: (i) inactivation of *NF2* in mouse
models led to a variety of malignant tumors but not to
mesothelioma [51], except if this inactivation was associ-
ated with asbestos exposure or with the inactivation of
other tumor suppressor genes [52]; recent studies
screening germline mutations in large cohorts of pa-
tients (reviewed in [44]) did not identify *NF2* as a cancer
susceptibility gene for MPM; and (iii) *NF2* mutations
showed a significantly higher mutation rate in MPM
with an advanced stage [22]. We previously showed that
MPM with mutations in members of the Hippo pathway
could be more sensitive to specific anti-cancer molecules
[53]. The Hippo pathway is becoming attractive for tar-
geted therapy in cancer, and numerous companies are
developing compounds to inhibit the Hippo pathway
[54]. Implementation of a therapeutic strategy based on
Hippo pathway dysregulation will need to take into ac-
count the clonality of *NF2* mutations in MPM. A hetero-
geneous frameshift deletion was also found in *CTNNB1*.
Inactivating mutations were previously reported in
MPM [22, 55]. Loss of the beta-catenin protein in a
MPM subclonal population lead to alterations in the
Wnt signaling pathway, which is involved in therapeutic
resistance [56]. Consequently, detection of this subclonal
population could be of therapeutic interest.

The presence of chromosomal abnormalities including
structural changes is a key feature of MPM [57]. How-
ever, our data and those of Chen et al. [18] showed that
chromosomal profiles are globally similar between differ-
ent regions within the same tumor. Similarly to copy
number aberrations, we did not identify specific relevant
gene fusions limited to one of the paired tumor samples.
As chromothripsis and chromoplexy lead to massive re-
arrangements in several chromosome regions in MPM
[58, 59], specific omic or NGS approaches will be needed
to deeply analyze spatial tumor heterogeneity at the
chromosome level. However, our results suggest that the
main chromosomal alterations occurring in MPM are
early events of mesothelial carcinogenesis, in accordance
with the mechanism of action of asbestos, which is well-
described as a genotoxic agent inducing chromosomal
damage [60].

In addition to the detection of heterogeneous popula-
tions characterized by specific mutations, our results
highlight major changes at the transcriptomic level be-
tween tumors at different anatomic sites in some pa-
tients, leading to the dysregulation of specific pathways.
We showed that gene expression changes are related to
epigenetic mechanisms in at least two patients. Methyla-
tion analysis is limited to 5 patients, with only two with
consistent changes, and allows only descriptive conclu-
sions to be drawn. Dysregulation of cell adhesion and
extracellular matrix organization pathways is frequent in

patients and may reflect variations between the propor-
tion of epithelioid and sarcomatoid phenotypes accord-
ing to the anatomical site. One histologic study showed
that a diagnosis of epithelioid MPM in the initial biopsy
was changed to the biphasic or sarcomatoid type in 19%
of cases when the surgical resection was evaluated [61].
In agreement, we observed variation greater than 10%
and up to 23% in the E/S.scores of the paired biopsies of
three patients [10]. We found epithelioid tumor samples
with high S.score (up to 0.68), consistent with previous
observations in larger series [10, 22]. This may be related
to the tumor representativity of a FFPE tumor section
compared to a piece of frozen tissue, or to a different
evaluation sensitivity. Another hypothesis is that the
S.score predicts the proportion of sarcomatoid-like cells,
which may not have completely a characteristic sarco-
matoid cell morphology. Gene expression-based signa-
tures to predict physiological process or prognosis such
as the E/S.score or the CV score (CLDN15/VIM) in
MPM [10, 11, 62], as well as response to treatment is be-
coming popular, and some of them have been evaluated
in phase 3 clinical trials for breast cancer [63]. Based on
our results, multi-site tumor sampling should be recom-
mended before implementing these assays in MPM.

Finally, we found consequent spatial intra-tumor het-
erogeneity of the immune microenvironment. We veri-
fied that this result was not due to a bias in sequencing
depth between samples. This is in line with previous
studies showing changes in the T cell repertoires at dif-
ferent anatomic sites [17, 18]. Using the MCP-counter
method, whose predictions were validated previously by
immunohistochemistry in MPM [10], we emphasized
differential immune cell infiltration between paired sam-
ples. Importantly, for three patients, the immune profile
could be considered hot or cold depending on the paired
samples. Substantial gene expression variation was also
observed for PDL1, PD1, and CTLA4 in some patients.
As immunotherapy with immune checkpoint inhibitors
is emerging as a promising therapeutic option for MPM
patients, biomarkers to predict the response to this
treatment are a crucial issue [2]. No predictive bio-
marker is currently clearly defined for MPM, but im-
mune cell infiltration and immune checkpoint
expression have been suggested in other cancers, indi-
cating that spatial intra-tumor heterogeneity needs to be
taken into account to identify biomarkers in MPM [64].

The limitation of this study is the small size of the
series in terms of patients and tumor samples per pa-
tient. Further studies in larger series are needed to
confirm the frequency of our major findings and to
have an exhaustive view of spatial intra-tumor hetero-
geneity in MPM. To predict response to therapy, it
would be also crucial to monitor clonal evolution
after treatment.
## Conclusions

Spatial intra-tumor heterogeneity is complex in MPM and varies among patients. We highlighted multiple types of heterogeneity, i.e., (i) genetic, (ii) transcriptomic, (iii) epigenetic, and (iv) linked to the immune micro-environment. It was found that the accuracy of histologic classification is increased by the examination of several tumor biopsies, and multi-sampling is recommended [21, 61]. Our molecular analysis also supports the notion that separate anatomic sites should be sampled from the pleural cavity to be able to estimate prognosis or predict response to treatment based on molecular characteristics, with the aim of developing molecular-based precision medicine strategies in order to improve patient survival and quality of life.

## Abbreviations

CCF: Cancer cell fraction; DE: Differentially expressed; DM: Differentially methylated; MPM: Malignant pleural mesothelioma; NGS: Next-generation sequencing; PET: Positron emission tomography; ssGSEA: Single-sample Gene Set Enrichment Analysis; TERT-prom: TERT promoter; WES: Whole exome sequencing

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s13073-021-00931-w.

**Additional file 1: Table S1.** Clinical annotations of MPM patients. **Table S2.** List of biomarkers included in MCP-counter. **Table S3.** Variants with structural consequences in genes and variants in the *TERT* promoter. **A** Germline mutations identified by whole exome sequencing. **B** Variants correspondence between sequencing methods. **C** Somatic variants with structural consequences in genes and variants in the *TERT* promoter identified by whole exome sequencing. **D** Variants with structural consequences in genes and variants in the *TERT* promoter identified by DNA targeted sequencing. **E** DNA targeted sequencing of multiple sections of the T227LE fresh tumor sample used to generate the MPM_83 primary cell line. **Table S4.** Clonality of the protein-altering somatic variants detected by whole exome sequencing. **A** Cancer cell fraction (CCF) values of protein-altering variants. **B** Clonality per patient based on protein-altering variants. **C** Cancer cell fraction (CCF) values of copy number variations (CNV) based on FACET analysis. **D** Clonality per patient based on FACET analysis. **Table S5.** Fusion transcripts identified by RNA-seq. **Table S6.** Pathway dysregulation between paired tumor samples. **A** Differentially expressed protein coding genes between paired tumor samples. **B** Pathways dysregulated between paired tumor samples in the GeneOntology Biological Process non-redundant database. **C** Pathways dysregulated between paired tumor samples in the KEGG database. **D** Pathways dysregulated between paired tumor samples in the Reactome database. **E** Pathways analysis using single sample Gene Set Enrichment Analysis (ssGSEA) in the GeneOntology Biological Process, KEGG and Reactome databases. **Table S7.** Epigenetic dysregulation between paired tumor samples. **A** Differentially methylated CpG between paired tumor samples. **B** Pathways dysregulated between paired tumor samples based on differentially expressed protein coding genes associated with differentially methylated CpG in the GeneOntology Biological Process non-redundant, KEGG and Reactome databases.

**Additional file 2: Figure S1.** Positron emission tomography (PET) scans of the two patients with biopsies performed on highest metabolite sites. **Figure S2.** Frequencies of *BAP1* and *NF2* variants in the fresh tumor sample used to generate the MPM_83 primary cell line. **Figure S3.** Cancer cell fraction in paired MPM samples. **Figure S4.** Pangenomic copy number variation profile determined by DNAcopy. **Figure S5.** Pangenomic copy number variation profile determined by FACET. **Figure S6.** Fusion transcripts identified by RNA-seq. **Figure S7.** Unsupervised clustering of MPM tumor samples based on gene expression profiles. **Figure S8.** Differentially expressed genes between paired tumor samples. **Figure S9.** Single sample Gene Set Enrichment Analysis (ssGSEA) for pathways related to cell adhesion and extracellular matrix. **Figure S10.** Epigenetic intra-tumor heterogeneity. **Figure S11.** Differential infiltration of stromal and immune cell populations. **Figure S12.** Distribution of immune and stromal populations among paired tumor samples. **Figure S13.** Single sample Gene Set Enrichment Analysis (ssGSEA) for pathways related to immune cell infiltration. **Figure S14.** Immunohistochemical staining of tumor-infiltrating T and B cells in patients with paired tumor samples switching from a cold to a hot immune profile.

## Acknowledgements

Not applicable.

## Authors' contributions

CM, JZR, FLPB, MCJ, and DJ are responsible for the study concept and design. MCC, LG, CB, EP, AS, and FLPB provided patient materials. JBA and FM collected and organized clinical data. MCC, LG, CB, and FGS performed a histologic analysis of tumor samples. SH and MCC performed immunohistochemistry analysis. LQ, JdW, JBA, FM, and CM contributed to sample preparation and quality control. CM performed the acquisition of data. YB, LM, QB, SC, TZH, and EL developed the bioinformatics tools and the pipeline analysis. SC, TZH, CM, and DJ performed the analysis and interpretation of data. CM, MCJ, and DJ were major contributors in writing the manuscript. FLPB, JZR, and DJ are responsible for the study supervision. All authors read and approved the final manuscript.

## Funding

This work was supported by Inserm, the Ligue Contre le Cancer (Ile de France committee and national program Cartes d'Identité des Tumeurs (CIT)), Hadassah France, the Chancellerie des Universités de Paris (Legs POIX) and Cancéropôle Ile-de-France (Emergence). JBA, FM, LQ and JdW were supported by grants from Fondation pour la Recherche Médicale (FRM), from the Institut thématique multi-organismes (ITMO) Cancer (Plan Cancer 2014-2019), from Cancéropôle Ile-de-France and from the Fondation ARC, respectively. JdW and FM were also supported by grant of the Société Française de Chirurgie Thoracique et Cardio-Vasculaire (SFCTCV). The Inserm research unit is supported by the Labex Oncolmmunology (investissement d'avenir), Coup d'Elan de la Fondation Bettencourt-Shueller, the SIRIC CARPEM, the Institut thématique multi-organismes (ITMO) Cancer (Plan Cancer 2014-2019) and Cancéropôle Ile-de-France.

## Availability of data and materials

The datasets supporting the conclusions of this article are included within the article (and its additional files). The NGS data (Targeted NGS, whole exome sequencing, RNA-Seq) are available in the European Genome-phenome Archive (EGA) under the EGAS00001005328 study number, https://ega.archive.org/studies/EGAS00001005328 [23]. The methylome data are available in Gene Expression Omnibus (GEO) repository under the GSE175769 series number, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE175769 [38].

## Declarations

### Ethics approval and consent to participate

This study is a part of a research project approved by a Medical Ethics Committee (CPP Ile-de-France II; ID-RCB: 2020-A01425-34). All patients gave their written informed consent for the use of their tumor specimen for research. The sampling procedures were approved by the French research ministry (CODECOH no. DC-2016-2771). The study methodologies were conformed to the standards set by the Declaration of Helsinki.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.
Author details

¹Centre de Recherche des Cordeliers, Inserm UMRS-1138, Sorbonne Université, Université de Paris, Functional Genomics of Solid Tumors, Paris, France. ²Present address: Service de Chirurgie Thoracique, Hôpital Calmette, CHRU de Lille, Lille, France. ³Present address: Service de Chirurgie Thoracique et Transplantation Pulmonaire, Hôpital Foch, Suresnes, France. ⁴University Paris-Est Créteil (UPEC), CEpiA (Clinical Epidemiology and Ageing), EA 7376-IMRB, UPEC, Créteil, France. ⁵GRC OncoThoParisEst, Service de Pneumologie, CHI Créteil, UPEC, Créteil, France. ⁶Programme Cartes d'identité des Tumeurs (CIT), Ligue Nationale Contre Le Cancer, Paris, France. ⁷Present address: IGDR UMR 6290, CNRS, Université de Rennes 1, Rennes, France. ⁸Assistance Publique-Hôpitaux de Paris, Hôpital Européen Georges Pompidou, Paris, France. ⁹Service d'Anatomopathologie et Cytologie, Université de Paris, Hôpital Européen Georges Pompidou, Paris, France. ¹⁰Univ. Lille, CHU Lille, Service de Chirurgie Thoracique, Hôpital Calmette, Lille, France. ¹¹Univ. Lille, CHU Lille, Institut de Pathologie, Lille, France. ¹²Univ. Lille, CNRS, Inserm, CHU Lille, Institut Pasteur de Lille, UMR9020 – UMR1277 – Canther – Cancer Heterogeneity, Plasticity and Resistance to Therapies, Lille, France. ¹³Centre National Référent MESOPATH, Centre Leon Berard, Lyon, France. ¹⁴Present address: Département de Pathologie Cellulaire et Tissulaire, CHU d'Angers, Angers, France. ¹⁵Univ. Lille, CHU Lille, Service de Pneumologie et d'Oncologie Thoracique, unité INSERM 1189 OncoThAl, Lille, France. ¹⁶Réseau National Expert pour le Méthéliome Pleural Malin (NETMESO), Lille, France. ¹⁷Service de Chirurgie Thoracique, Hôpital Européen Georges Pompidou, Paris, France.

Received: 9 November 2020 Accepted: 30 June 2021

Published online: 14 July 2021

References

1. Zalcman G, Mazieres J, Margery J, Greillier L, Audigier-Valette C, Moro-Sibilot D, et al. Bevacizumab for newly diagnosed pleural mesothelioma in the Mesothelioma Avastin Cisplatin Pemetrexed Study (MAPS): a randomised, controlled, open-label, phase 3 trial. Lancet. 2016;387(10026):1405-14. https://doi.org/10.1016/S0140-6736(15)01238-6.

2. de Gooijer CJ, Borm FJ, Scherpereel A, Baas P. Immunotherapy in malignant pleural mesothelioma. Front Oncol. 2020;10:187. https://doi.org/10.3389/fonc.2020.00187.

3. Scherpereel A, Mazieres J, Greillier L, Lantuejoul S, Dó P, Bylicki O, et al. Nivolumab or nivolumab plus ipilimumab in patients with relapsed malignant pleural mesothelioma (IFCT-1501 MAPS2): a multicentre, open-label, randomised, non-comparative, phase 2 trial. Lancet Oncol. 2019;20(2): 239-53. https://doi.org/10.1016/S1470-2045(18)30765-4.

4. Baas P, Scherpereel A, Nowak AK, Fujimoto N, Peters S, Tsao AS, et al. First-line nivolumab plus ipilimumab in unresectable malignant pleural mesothelioma (CheckMate 743): a multicentre, randomised, open-label, phase 3 trial. The Lancet. 2021;397(10272):375-86. https://doi.org/10.1016/S0140-6736(20)32714-8.

5. Oehl K, Vrugt B, Opitz I, Meerang M. Heterogeneity in malignant pleural mesothelioma. Int J Mol Sci. 2018;19(6):1603. https://doi.org/10.3390/ijms19061603.

6. Husain AN, Colby TV, Ordóñez NG, Allen TC, Attanoos RL, Beasley MB, et al. Guidelines for Pathologic Diagnosis of Malignant Mesothelioma 2017 Update of the Consensus Statement From the International Mesothelioma Interest Group. Arch Pathol Lab Med. 2018;142(1):89-108. https://doi.org/10.5858/arpa.2017-0124-RA.

7. Husain AN, Colby T, Ordonez N, Krausz T, Attanoos R, Beasley MB, et al. Guidelines for pathologic diagnosis of malignant mesothelioma: 2012 update of the consensus statement from the International Mesothelioma Interest Group. Arch Pathol Lab Med. 2013;137(5):647-67. https://doi.org/10.5858/arpa.2012-0214-OA.

8. Bilecz A, Stockhammer P, Theegarten D, Kern I, Jakopovic M, Samarzija M, et al. Comparative analysis of prognostic histopathologic parameters in subtypes of epithelioid pleural mesothelioma. Histopathology. 2020;77(1): 55-66. https://doi.org/10.1111/his.14105.

9. Galateau Salle F, Le Stang N, Tirode F, Courtiol P, Nicholson AG, Tsao M-S, et al. Comprehensive molecular and pathologic evaluation of transitional mesothelioma assisted by deep learning approach: a multi-institutional study of the international mesothelioma panel from the MESOPATH reference center. J Thorac Oncol. 2020;15(6):1037-53. https://doi.org/10.1016/j.jtho.2020.01.025.

10. Blum Y, Meiller C, Quetel L, Elarouci N, Ayadi M, Tashtanbaeva D, et al. Dissecting heterogeneity in malignant pleural mesothelioma through histomolecular gradients for clinical applications. Nat Commun. 2019;10(1):1333. https://doi.org/10.1038/s41467-019-09307-6.

11. Bueno R, Stawiski EW, Goldstein LD, Durinck S, De Rienzo A, Modrusan Z, et al. Comprehensive genomic analysis of malignant pleural mesothelioma identifies recurrent mutations, gene fusions and splicing alterations. Nat Genet. 2016;48(4):407-16. https://doi.org/10.1038/ng.3520.

12. de Reyniès A, Jaurand M-C, Renier A, Couchy G, Hysi I, Elarouci N, et al. Molecular classification of malignant pleural mesothelioma: identification of a poor prognosis subgroup linked to the epithelial-to-mesenchymal transition. Clin Cancer Res. 2014;20(5):1323-34. https://doi.org/10.1158/1078-0432.CCR-13-2429.

13. Hmeljak J, Sanchez-Vega F, Hoadley KA, Shih J, Stewart C, Heiman D, et al. Integrative molecular characterization of malignant pleural mesothelioma. Cancer Discov. 2018;8(12):1548-65. https://doi.org/10.1158/2159-8290.CD-18-0804.

14. Alcala N, Mangiante L, Le-Stang N, Gustafson CE, Boyault S, Damiola F, et al. Redefining malignant pleural mesothelioma types as a continuum uncovers immune-vascular interactions. EBioMedicine. 2019;48:191-202. https://doi.org/10.1016/j.ebiom.2019.09.003.

15. Comertpay S, Pastorino S, Tanji M, Mezzapelle R, Christiane O, Napolitano A, et al. Evaluation of clonal origin of malignant mesothelioma. J Transl Med. 2014;12(1):301. https://doi.org/10.1186/s12967-014-0301-3.

16. Oey H, Daniels M, Relan V, Chee TM, Davidson MR, Yang IA, et al. Whole-genome sequencing of human malignant mesothelioma tumours and cell lines. Carcinogenesis. 2019;40(6):724-34. https://doi.org/10.1093/carcin/bgz066.

17. Kiyotani K, Park J-H, Inoue H, Husain A, Olugbile S, Zewde M, et al. Integrated analysis of somatic mutations and immune microenvironment in malignant pleural mesothelioma. Oncoimmunology. 2017;6(2):e1278330. https://doi.org/10.1080/2162402X.2016.1278330.

18. Chen R, Lee W-C, Fujimoto J, Li J, Hu X, Mehran R, et al. Evolution of genomic and T-cell repertoire heterogeneity of malignant pleural mesothelioma under dasatinib treatment. Clin Cancer Res. 2020;26(20): 5477-86. https://doi.org/10.1158/1078-0432.CCR-20-1767.

19. Minnema-Luiting J, Vroman H, Aerts J, Cornelissen R. Heterogeneity in immune cell content in malignant pleural mesothelioma. Int J Mol Sci. 2018; 19(4). https://doi.org/10.3390/ijms19041041.

20. Galateau-Sallé F, Gilg Soit Ilg A, Le Stang N, Brochard P, Pairon JC, Astoul P, et al. The French mesothelioma network from 1998 to 2013. Ann Pathol. 2014;34(1):51-63. https://doi.org/10.1016/j.annpat.2014.01.009.

21. Nicholson AG, Sauter JL, Nowak AK, Kindler HL, Gill RR, Remy-Jardin M, et al. EURACAN/IASLC Proposals for updating the histologic classification of pleural mesothelioma: towards a More multidisciplinary approach. J Thorac Oncol. 2020;15(1):29-49. https://doi.org/10.1016/j.jtho.2019.08.2506.

22. Quetel L, Meiller C, Assié J-B, Blum Y, Imbeaud S, Montagne F, et al. Genetic alterations of malignant pleural mesothelioma: association to tumor heterogeneity and overall survival. Mol Oncol. 2020;14(6):1207-23. https://doi.org/10.1002/1878-0261.12651.

23. EGA study EGAS00001005328. https://ega-archive.org/studies/EGAS00001005328.

24. Haeussler M, Zweig AS, Tyner C, Speir ML, Rosenbloom KR, Raney BJ, et al. The UCSC Genome Browser database: 2019 update. Nucleic Acids Res. 2019; 47(D1):D853-8. https://doi.org/10.1093/nar/gky1095.

25. Cibulskis K, Lawrence MS, Carter SL, Sivachenko A, Jaffe D, Sougnez C, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat Biotechnol. 2013;31(3):213-9. https://doi.org/10.1038/nbt.2514.

26. Tarasov A, Vilella AJ, Cuppen E, Nijman UJ, Prins S. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31(12):2032-4. https://doi.org/10.1093/bioinformatics/btv098.

27. Venkatraman ES, Olshen AB. A faster circular binary segmentation algorithm for the analysis of array CGH data. Bioinformatics. 2007;23(6):657-63. https://doi.org/10.1093/bioinformatics/btl646.

28. Shen R, Seshan VE. FACETS: allele-specific copy number and clonal heterogeneity analysis tool for high-throughput DNA sequencing. Nucleic Acids Res. 2016;44(16):e131. https://doi.org/10.1093/nar/gkw520.

29. Sondka Z, Bamford S, Cole CG, Ward SA, Dunham I, Forbes SA. The COSMIC Cancer Gene Census: describing genetic dysfunction across all human
cancers. Nat Rev Cancer. 2018;18(11):696-705. https://doi.org/10.1038/s41
568-018-0060-1.

30. Cosmic Cancer Gene Census. http://cancer.sanger.ac.uk/census. Accessed 24 Jan 2020.

31. Shinde J, Bayard Q, Imbeaud S, Hirsch TZ, Liu F, Renault V, et al. Palimpsest: an R package for studying mutational and structural variant signatures along clonal evolution in cancer. Bioinformatics. 2018;34(19):3380-1. https://doi.org/10.1093/bioinformatics/bty388.

32. Letouzé E, Shinde J, Renault V, Couchy G, Blanc J-F, Tubacher E, et al. Mutational signatures reveal the dynamic interplay of risk factors and cellular processes during liver tumorigenesis. Nat Commun. 2017;8(1):1315. https://doi.org/10.1038/s41467-017-01358-x.

33. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36. https://doi.org/10.1186/gb-2013-14-4-r36.

34. Anders S, Pyl PT, Huber W. HTSeq-a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166-9. https://doi.org/10.1093/bioinformatics/btu638.

35. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. https://doi.org/10.1186/s13059-014-0550-8.

36. Sun L, Dong S, Ge Y, Fonseca JP, Robinson ZT, Mysore KS, et al. DiVenn: An Interactive and Integrated Web-Based Visualization Tool for Comparing Gene Lists. Front Genet. 2019;10:421. https://doi.org/10.3389/fgene.2019.00421.

37. Bayard Q, Meunier L, Peneau C, Renault V, Shinde J, Nault J-C, et al. Cyclin A2/E1 activation defines a hepatocellular carcinoma subclass with a rearrangement signature of replication stress. Nat Commun. 2018;9(1):5235. https://doi.org/10.1038/s41467-018-07552-9.

38. GEO series GSE175769. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?aacc=GSE175769.

39. Liao Y, Wang J, Jaehnig EJ, Shi Z, Zhang B. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res. 2019;47(W1): W199-205. https://doi.org/10.1093/nar/gkz401.

40. Häzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14(1):7. https://doi.org/10.1186/1471-2105-14-7.

41. Becht V, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016;17(1): 218. https://doi.org/10.1186/s13059-016-1070-5.

42. EGA study EGAS00001001563. https://ega-archive.org/studies/EGAS00001001563.

43. Broad Institute TCGA Genome Data Analysis Center. Analysis-ready standardized TCGA data from Broad GDAC Firehose stddata__2015_06_01 run. 2015. https://doi.org/10.7908/C1251HBG.

44. Panou V, Røe OD. Inherited Genetic Mutations and Polymorphisms in Malignant Mesothelioma: A Comprehensive Review. JMS. 2020;21(12):4327. https://doi.org/10.3390/ijms21124327.

45. Chen Y, Tseng S-H. Targeting tropomyosin-receptor kinase fused gene in cancer. Anticancer Res. 2014;34(4):1595-600.

46. Birkbak NJ, McGranahan N. Cancer genome evolutionary trajectories in metastasis. Cancer Cell. 2020;37(1):8-19. https://doi.org/10.1016/j.ccell.2020.09.12004.

47. Amintas S, Bedel A, Moreau-Gaudry F, Boutin J, Buscail L, Merlio J-P, et al. Circulating tumor cell clusters: united we stand divided we fall. JIMS. 2020; 21(7):2653. https://doi.org/10.3390/ijms21072653.

48. Hjerpe A, Ascoli V, Bedrossian CWM, Boon ME, Creaney J, Davidson B, et al. Guidelines for the cytopathologic diagnosis of epithelioid and mixed-type malignant mesothelioma: a secondary publication. Cytopathology. 2015; 26(3):142-56. https://doi.org/10.1111/cyt.12250.

49. Matsumoto S, Nabeshima K, Kamei T, Hiroshima K, Kawahara K, Hata S, et al. Morphology of 9p21 homozygous deletion-positive pleural mesothelioma cells analyzed using fluorescence in situ hybridization and virtual microscope system in effusion cytology: 9p21 Deletion-Positive MPM Cells. Cancer Cytopathol. 2013;121(8):415-22. https://doi.org/10.1002/cncy.21269.

50. Sato T, Sekido Y. NF2/Merlin inactivation and potential therapeutic targets in mesothelioma. JIMS. 2018;19(4):988. https://doi.org/10.3390/ijms19040988.

51. McClatchey AI, Saotome I, Mercer K, Crowley D, Gusella JF, Bronson RT, et al. Mice heterozygous for a mutation at the NF2 tumor suppressor locus

develop a range of highly metastatic tumors. Genes Dev. 1998;12(8):1121-33. https://doi.org/10.1101/gad.12.8.1121.

52. Jean D, Jaurand M-C. Mesotheliomas in Genetically Engineered Mice Unravel Mechanism of Mesothelial Carcinogenesis. JMS. 2018;19(8):2191. https://doi.org/10.3390/ijms19082191.

53. Tranchant R, Quetel L, Tallet A, Meiller C, Renier A, de Koning L, et al. Co-occurring mutations of tumor suppressor genes, LATS2 and NF2, in malignant pleural mesothelioma. Clin Cancer Res. 2017;23(12):3191-202. https://doi.org/10.1158/1078-0432.CCR-16-1971.

54. Dey A, Varelas X, Guan K-L. Targeting the Hippo pathway in cancer, fibrosis, wound healing and regenerative medicine. Nat Rev Drug Discov. 2020;19(7): 480-94. https://doi.org/10.1038/s41573-020-0070-z.

55. Shigemitsu K, Sekido Y, Usami N, Mori S, Sato M, Horio Y, et al. Genetic alteration of the beta-catenin gene (CTNNB1) in human lung cancer and malignant mesothelioma and identification of a new 3p21.3 homozygous deletion. Oncogene. 2001;20(31):4249-57. https://doi.org/10.1038/sj.onc.1204557.

56. Martin-Orozco E, Sanchez-Fernandez A, Ortiz-Parra I, Ayala-San Nicolas M. WNT signaling in tumors: the way to evade drugs and immunity. Front Immunol. 2019;10:2854. https://doi.org/10.3389/fimmu.2019.02854.

57. Blanquart C, Jaurand M-C, Jean D. The biology of malignant mesothelioma and the relevance of preclinical models. Front Oncol. 2020;10:388. https://doi.org/10.3389/fonc.2020.00388.

58. Mansfield AS, Peikert T, Smadbeck JB, Udell JBM, Garcia-Rivera E, Elsbernd L, et al. Neoantigenic potential of complex chromosomal rearrangements in mesothelioma. J Thorac Oncol. 2019;14(2):276-87. https://doi.org/10.1016/j.jtho.2018.10.001.

59. Yoshikawa Y, Emi M, Hashimoto-Tamaoki T, Ohmuraya M, Sato A, Tsujimura T, et al. High-density array-CGH with targeted NGS unmask multiple noncontiguous minute deletions on chromosome 3p21 in mesothelioma. Proc Natl Acad Sci USA. 2016;113(47):13432-7. https://doi.org/10.1073/pnas.1612074113.

60. Kane AB, Jean D, Knuutila S, Jaurand M-C. Malignant Mesothelioma: Mechanism of Carcinogenesis. In: Anttila S, Boffetta P, editors. Occupational Cancers [Internet]. Cham: Springer International Publishing; 2020. p. 343-62. [cited 2020 Feb 15]. Available from: http://link.springer.com/10.1007/978-3-030-30766-0_19.

61. Chirieac LR, Hung YP, Foo WC, Hofer MD, VanderLaan PA, Richards WG, et al. Diagnostic value of biopsy sampling in predicting histology in patients with diffuse malignant pleural mesothelioma. Cancer. 2019;125(23):4164-71. https://doi.org/10.1002/cncr.32416.

62. Severson DT, De Rienzo A, Bueno R. Mesothelioma in the age of "Omics": Before and after The Cancer Genome Atlas. J Thorac Cardiovasc Surg. 2020; 160:1078-1083.e2.

63. Markopoulos C, Hyams DM, Gomez HL, Harries M, Nakamura S, Traina T, et al. Multigene assays in early breast cancer: Insights from recent phase 3 studies. Eur J Surg Oncol. 2020;46(4):656-66. https://doi.org/10.1016/j.ejso.2019.010.19.

64. Bai R, Lv Z, Xu D, Cui J. Predictive biomarkers for cancer immunotherapy with immune checkpoint inhibitors. Biomark Res. 2020;8(1):34. https://doi.org/10.1186/s40364-020-00209-0.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_9.png" /></figure>