Redefining malignant pleural mesothelioma types as a continuum uncovers immune-vascular interactions

Abstract

Background: Malignant Pleural Mesothelioma (MPM) is an aggressive disease related to asbestos exposure, with no effective therapeutic options.

Methods: We undertook unsupervised analyses of RNA-sequencing data of 284 MPMs, with no assumption of discreteness. Using immunohistochemistry, we performed an orthogonal validation on a subset of 103 samples and a biological replication in an independent series of 77 samples.

Findings: A continuum of molecular profiles explained the prognosis of the disease better than any discrete model. The immune and vascular pathways were the major sources of molecular variation, with strong differences in the expression of immune checkpoints and pro-angiogenic genes; the extrema of this continuum had specific molecular profiles: a "hot" bad-prognosis profile, with high lymphocyte infiltration and high expression of immune checkpoints and pro-angiogenic genes; a "cold" bad-prognosis profile, with low lymphocyte infiltration and high expression of pro-angiogenic genes; and a "VEGFR2+/VISTA+" better-prognosis profile, with high expression of immune checkpoint VISTA and pro-angiogenic gene VEGFR2. We validated the gene expression levels at the protein level for a subset of five selected genes belonging to the immune and vascular pathways (CD8A, PDL1, VEGFR3, VEGFR2, and VISTA), in the validation series, and replicated the molecular profiles as well as their prognostic value in the replication series.

Interpretation: The prognosis of MPM is best explained by a continuous model, which extremes show specific expression patterns of genes involved in angiogenesis and immune response.

Keywords: Angiogenesis; French MESOBANK; Immunotherapy; MESOMICS project; Pleural mesothelioma.

Fig. 1

Malignant Pleural Mesothelioma expression profiles follow a continuum model. a) Two-dimensional summary of 284 transcriptomes using Principal Component Analysis (PCA). Point colours represent the three histopathological types, and the overlayed blue-coloured rectangles represent the survival in nine regions; the filled shapes on the bottom panel correspond to the density of samples from each histopathological type on Dimension 1, and the filled shapes on the right panel correspond to the RNA-seq-estimated mean proportion of immune cells from 10 cell types, in each sample, as a function of Dimension 2 coordinates, computed using a moving average with a window size of 30 Dimension 2 units. b) Integral AUC (iAUC) of five Cox proportional hazards survival models: (i) a model based on the three histological types; (ii) a model based on the percentage of sarcomatoid; (iii) a model based on the four molecular clusters from the study of Bueno and colleagues ; (iv) a model based on the coordinates of samples on Dimension 1; (v) a model based on the coordinates of samples on Dimensions 1 and 2. c) Gene-Set Enrichment Analysis (GSEA) of the genes most correlated with Dimensions 1 and 2, based on the hallmarks of cancer gene sets; violin plots and boxplots represent the distribution of Pearson correlation coefficients between gene expression and Dimensions 1 (red) or 2 (blue); genes in parenthesis are not part of the current hallmark annotation; only the three hallmarks with the highest correlation are represented; see Fig. S9 for the results of all hallmarks. In the boxplot representation, centre line represents the median and box bounds represent the inter-quartile range (IQR). The whiskers span a 1.5-fold IQR or the highest and lowest observation values if they extend no further than the 1.5-fold IQR. d) Correlation circle of the Principal Component Analysis (PCA) from panel (a) for 12 genes of interest. Arrow lengths and direction correspond to the strength and sign of the correlation between the variable and Dimensions 1 and 2. e) Forest plot of hazard ratios for overall survival with age, sex and asbestos exposure as covariables. The black boxes represent estimated hazard ratios and whiskers represent the associated 95% confidence intervals. Wald test q-values are shown on the right. Only the markers significantly associated with survival are represented (Wald test q < 0.05); see Table S8 for the results of all genes. Data used in (a) and (c) correspond to the n = 211 samples from the study of Bueno and colleagues and the n = 73 transcriptomes from the TCGA MESO cohort . Data used in (b) correspond to the n = 199 samples from the Bueno cohort with RNA-seq data and available percentage of sarcomatoid component. Data used in (e) correspond to n = 205 samples from the Bueno cohort and n = 59 samples from the TCGA MESO cohort with RNA-seq data and available asbestos exposure annotations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

Technical validation of a five-gene panel on 103 MPMs. a) Left panel: correlation circle of the Principal Component Analysis (PCA) based on the RNA-seq expression of the five-gene panel. Arrow lengths and direction correspond to the strength and sign of the correlation between the variable and Dimensions 1 and 2. Data used correspond to the n = 211 samples from the study of Bueno and colleagues and the n = 73 transcriptomes from the TCGA MESO cohort . Right panel: correlation matrix of the five-gene panel expression (upper triangle), of their protein expression (lower triangle), and correlation between expression from RNA-seq data and protein expression from IHC data (green diagonal). Colours correspond to the magnitude and sign of the correlations and statistically significant correlations are surrounded by a black box; dendrograms represent hierarchical clustering of gene or protein expression levels. Data used correspond to the n = 103 samples from the Bueno cohort with RNA-seq data and with Tissue MicroArray (TMA) IHC staining data. b) Tissue MicroArray (TMA) IHC staining from the technical validation series corresponding to n = 103 samples from the Bueno cohort , with 0.6 mm core diameter at 5.2× magnification, for the five-gene panel, representing the positive and negative references of the tested protein expression. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

Diagnostic and prognostic value of the expression level of the five-protein panel in the replication series of 77 MPMs, determined by immunohistochemistry. a) Top panel: two-dimensional summary of the gene expression in the discovery cohort (n = 284) (PCA; subset of Fig. 1a). Bottom panel: two-dimensional summary of the protein expression of the five genes in the replication cohort (n = 77) (IHC PCA). Point colours correspond to the three sample sets from Table 1. b) Top panel: correlation circle of the IHC PCA (n = 77) from (a) bottom panel, where arrow lengths and direction correspond to the strength and sign of the correlation between the variable and Dimensions 1 and 2. Bottom panel: correlation matrix of the protein expression of the 77 MPMs from the replication cohort, where colours correspond to the magnitude and sign of the correlations; the dendrogram represents a hierarchical clustering of protein expression. Significant correlations are surrounded by a black box. c) Left panel: gene expression levels (normalized read counts) in the discovery cohort between long-survival epithelioid, long-survival epithelioid, and sarcomatoid groups, resulting in n = 82 samples from the study of Bueno and colleagues and the n = 31 transcriptomes from the TCGA MESO cohort , for the three sets; each row presents violin plots and boxplots for a gene, with stars representing the significance level of pairwise comparisons between groups (q-values from two-sided independent Wilcoxon U tests). Right panel: Protein expression levels (% of cells where the protein is expressed) in the replication cohort, for the three sets; each row presents violin plots and boxplots for a protein, with stars representing the significance level of pairwise comparisons between groups (q-values from one-sided paired Wilcoxon T-tests). Sample sizes were 74, 74, 73, 70, 74, and 69, for CD8, PDL1 in the tumour, PDL1 in TILS, VEGFR2, VEGFR3, and VISTA, respectively. In the boxplot representation, centre line represents the median and box bounds represent the inter-quartile range (IQR). The whiskers span a 1.5-fold IQR or the highest and lowest observation values if they extend no further than the 1.5-fold IQR. d) PDL1 immunohistochemistry of two MPM cases from the replication cohort (left panel: short-survival epithelioid sample; right panel: sarcomatoid sample), both PDL1+ and PDL1 TILS+. Upper panels: Hematoxylin Eosin Saffron (HE) stain at 7× magnification, where white and black arrows show tumour cells and TILS, respectively. Lower panels: corresponding staining with PDL1 rabbit monoclonal antibody (cl SP263; VENTANA) at 7× magnification, where white and black arrows show positive staining of tumour cells and TILS, respectively. e) Protein expression level of VISTA in the replication cohort when considering epithelioid subtypes, independently of the sample set (upper panel) and in addition to the sample set (bottom panel). Data used correspond to n = 63 samples from the replication cohort with available data for all protein markers. In the boxplot representation, centre line represents the median and box bounds represent the inter-quartile range (IQR). The whiskers span a 1.5-fold IQR or the highest and lowest observation values if they extend no further than the 1.5-fold IQR.

Fig. 4

Characteristics of the three Malignant Pleural Mesothelioma transcriptomic profiles. The schematic position of samples harbouring a given profile in the two-dimensional summary from Fig. 1a (n = 284) is represented in the bottom right panel. For each profile, the hallmarks of cancer generally upregulated are indicated by pictograms in the upper left part, the histological type composition is represented by a pie chart in the upper right part, the proportion of tumour infiltrating lymphocytes estimated from the RNA-seq data (Figs. 1 and S10) is represented by a bar plot in the bottom left part, and the expression of representative genes is represented by a radar plot in the bottom right part. Tissue MicroArray (TMA) IHC staining from the technical validation series, with 0.6 mm core diameter at 5.2× magnification, for the five-gene panel, are presented above each panel.