Analysis of PD-L1 and PD-L2 expression correlated genes in CCLE dataset
We first addressed PD-L1 and PD-L2 mRNA expression in 114 NSCLC cell lines from the CCLE project Lung_NSC. Retrospective data analyses showed that PD-L1 and PD-L2 mRNA expression was correlated (Fig. 1b, c). PD-L2 mRNA expression was in general 2- to 4-fold lower than PD-L1 mRNA expression (Fig. 1b). Neither PD-L1 nor PD-L2 expression correlated with their respective receptors PD-1, CD80, and RGMB with the criteria for significant expression correlation being a Spearman correlation coefficient r ≥ 0.4 or ≤ − 0.4, a Pearson correlation coefficient r ≥ 0.3 or ≤ − 0.3, and all P values < 0.05 (Fig. 1c). Notably, correlation between RGMB and PD-L1 expression, as well as between RGMB and PD-L2 expression, was near significance (Fig. 1c). We next identified PD-L1 and PD-L2 expression correlated genes in the CCLE dataset Lung_NSC. 489 genes expression correlated with PD-L1, 191 genes expression correlated with PD-L2, and 111 genes expression correlated with both PD-L1 and PD-L2 (Fig. 1d and Additional file 1: Table S1). Interestingly, GSEA revealed that PD-L1 and PD-L2 expression correlated genes were enriched to the hallmark gene sets TNFα signaling via Nuclear factor kappa B (NFKB) (genes regulated by NF-kB in response to TNF), Kirsten rat sarcoma viral oncogene homolog (KRAS) signaling (genes up-regulated by KRAS activation), and Transforming growth factor β (TGFβ) signaling (genes up-regulated in response to TGFB1) (Additional file 2: Fig. S1B and Additional file 3: Table S2). GSEA also showed that PD-L1 and PD-L2 expression correlated genes were enriched for belonging to the epithelial-mesenchymal-transition (EMT) hallmark gene set (genes defining the epithelial-mesenchymal-transition, as in wound healing, fibrosis, and metastasis) (Additional file 2: Fig. S1B and Additional file 3: Table S2). We note a previously described association between EMT, PD-L1 expression, and immunosuppression [38, 39]. GSEA showed that among the hallmark gene sets enriched with PD-L1 expression correlated genes, but not with PD-L2 expression correlated genes, were the hallmark gene sets IFN-α response (genes up-regulated in response to alpha interferon proteins) and IFN-γ response (genes up-regulated in response to IFNγ) (Additional file 3: Table S2). We next addressed how individual genes involved in IFN signaling expression correlated with PD-L1 and PD-L2 in Lung_NSC (Additional file 4: Table S3). Genes in the canonical signaling pathway for IFN-γ include IFN-γ receptor 1 and 2 (IFNGR1 and IFNGR2), Janus kinase 1 and 2 (JAK1 and JAK2), Signal transducer and activator of transcription 1, 2, and 3 (STAT1, STAT2, and STAT3), and Interferon regulatory factor 1 (IRF1). Genes in the canonical signaling pathway for IFN-α include IFN-α receptor 1 and 2 (IFNAR1 and IFNAR2), JAK1 and Tyrosine kinase 2 (TYK2), STAT1, STAT2, and STAT3, and Interferon regulatory factor 9 (IRF9) [21]. In Lung_NSC, 152 genes expression correlated with IRF1 and 248 genes expression correlated with IRF9 (Additional file 2: Fig. S1A and Additional file 5: Table S4). Interestingly, IRF9, but not IRF1, expression correlated with PD-L1 (Additional file 4: Table S3). IRF1 and IRF9 expression correlated not with PD-L2 (Additional file 4: Table S3). We also examined whether the identified Lung_NSC PD-L1 and PD-L2 expression correlated genes (Additional file 1: Table S1) were represented among the identified Lung-NSC IRF1 and IRF9 expression correlated genes (Additional file 5: Table S4). Of the PD-L1 expression correlated genes, 10% (49/489) expression correlated with IRF1, and 17% (83/489) expression correlated with IRF9 (Additional file 2: Fig. S1A). Of the PD-L2 expression correlated genes, 6% (11/191) expression correlated with IRF1, and 5% (9/191) expression correlated with IRF9 (Additional file 2: Fig. S1A). IRF1 and IRF9 were not expression correlated with PD-1, CD80, and RGMB (Additional file 4: Table S3). GSEA showed extensive overlaps between hallmark gene sets enriched with PD-L1 and PD-L2 expression correlated genes and hallmark gene sets enriched with IRF1 and IRF9 expression correlated genes (Additional file 2: Fig. S1B, Additional file 3: Table S2, and Additional file 6: Table S5).
PD-L1 and PD-L2 expression correlated genes are differently associated with IFN signaling and immune cell markers in NSCLC tumor samples
To investigate PD-L1 and PD-L2 mRNA expression in NSCLC tumor samples we analyzed available RNA sequencing data from the TCGA projects LUAD and LUSC datasets, as well as RNA sequencing data for normal lung tissue from the GTEx project. In the LUAD and LUSC datasets, PD-L1 expression was 2- to 3-fold higher than PD-L2 expression (Additional file 7: Table S6). PD-L1 and PD-L2 expression correlated in LUSC and LUAD tumor samples, and as well as in normal lung tissue (Fig. 2a–c). In LUAD, both PD-1 and CD80 expression correlated with PD-L1, as well as with PD-L2 (Fig. 2b). In LUSC, CD80 expression correlated with PD-L2 (Fig. 2c). In the LUAD and LUSC datasets, no expression correlation between RGMB and PD-L1 or between RGMB and PD-L2 was identified (Fig. 2b, c). We note that CD80 expression correlated with PD-1 in both the LUAD and LUSC datasets (Fig. 2b, c).
Next, we identified PD-L1 and PD-L2 expression correlated genes in LUAD and LUSC datasets. For the LUAD dataset, 257 genes and 914 genes expression correlated with PD-L1 and PD-L2, respectively (Fig. 2d and Additional file 1: Table S1). We noted that 211 genes expression correlated with both PD-L1 and PD-L2, representing 82% of the PD-L1 expression correlated genes and 23% of the PD-L2 expression correlated genes (Fig. 2d). GSEA showed that PD-L1 and PD-L2 expression correlated genes were enriched for belonging to the same hallmark gene sets, except for the EMT gene set only being significant among the PD-L2 expression correlated genes (Additional file 3: Table S2 and Additional file 8: Fig. S2B). Notably, hallmark gene sets IFN-α response and IFN-γ response were among the hallmark gene sets for both PD-L1 and PD-L2 expression correlated genes (Additional file 3: Table S2 and Additional file 8: Fig. S2B). In the LUSC dataset, 26 genes expression correlated with PD-L1, and 326 genes expression correlated with PD-L2 (Fig. 2e and Additional file 1: Table S1). Of these, 13 genes expression correlated with both PD-L1 and PD-L2, representing 50% of the PD-L1 expression correlated genes but only 4% of the PD-L2 expression correlated genes (Fig. 2e). GSEA showed that all hallmark gene sets enriched for PD-L1 expression correlated genes also were enriched for PD-L2 expression correlated genes (Additional file 3: Table S2 and Additional file 8: Fig. S2D).
We also investigated the connection between PD-L1 and PD-L2 expression correlated genes and IFN signaling, with IFN signaling being represented by the 13 genes IFNG (encoding IFN-γ), IFNGR1, IFNGR2, IFNAR1, IFNAR2, JAK1, JAK2, TYK2, STAT1, STAT2, STAT3, IRF1, and IRF9 (Additional file 9: Table S7). We performed an expression correlation analysis between these 13 genes and PD-L1 and PD-L2 (Additional file 4: Table S3). In the LUAD dataset, PD-L1 and PD-L2 expression correlated with IRF1, IFNG, STAT1, and JAK2 (Additional file 4: Table S3). In addition, PD-L2 expression correlated with IFNGR1 (Additional file 4: Table S3). Notably, PD-L1 expression was not correlated with IRF9, a different scenario from that observed in CCLE Lung_NSC. In line with this, we noted that in the LUAD dataset the number of expression correlated genes with IRF1 (559 genes) was higher than that for IFR9 (82 genes) (Additional file 5: Table S4 and Additional file 8: Fig. S2A) which differed from the observations for CCLE Lung_NSC (Additional file 2: Fig. S1A and Additional file 5: Table S4). Prominent numbers of commonly correlated genes for IRF1 and PD-L1 (n = 141), as well as for IRF1 and PD-L2 (n = 488), were observed in LUAD (Additional file 8: Fig. S2A). This further illustrated that only 66/559 IRF1 expression correlated genes were not PD-L1 or PD-L2 expression correlated (Additional file 8: Fig. S2A). In comparison, there were fewer commonly correlated genes for IRF9 and PD-L1 (n = 19) and for IRF9 and PD-L2 (n = 44) (Additional file 8: Fig. S2A). GSEA showed extensive overlapping among hallmark gene sets enriched with PD-L1, PD-L2, IRF1, and IRF9 expression correlated genes (Additional file 3: Table S2, Additional file 6: Table S5 and Additional file 8: Fig. S2B).
In the LUSC dataset, of the 13 IFN signaling genes JAK2 and IRF1 expression correlated with PD-L1 (Additional file 4: Table S3). Moreover, 10/26 of the PD-L1 expression correlated genes were IRF1 expression correlated and 3/26 of the PD-L1 expression correlated genes were IRF9 expression correlated (Additional file 8: Fig. S2C). GSEA showed that the hallmark gene sets IFN-γ response and allograft rejection (genes up-regulated during transplant rejection) were common hallmark gene sets for both PD-L1 and IRF1 expression correlated genes, as well as for PD-L1 and IRF9 expression correlated genes (Additional file 3: Table S2, Additional file 6: Table S5 and Additional file 8: Fig. S2D). In the LUSC dataset, of the 13 IFN signaling genes IFNG, JAK2, STAT1, and IRF1 expression correlated with PD-L2 (Additional file 4: Table S3). As also observed in the LUAD dataset, in the LUSC dataset the number of IRF1 expression correlated genes (n = 751) was more pronounced than the number of IFR9 expression correlated genes (n = 90) (Additional file 5: Table S4 and Additional file 8: Fig. S2A, C). 297/326 of the PD-L2 expression correlated genes were IRF1 expression correlated and 18/326 of the PD-L2 expression correlated genes were IRF9 expression correlated (Additional file 8: Fig. S2C). GSEA showed that PD-L2, IRF1, and IRF9 expression correlated genes to a large extent belong to the same hallmark gene sets (Additional file 3: Table S2, Additional file 6: Table S5 and Additional file 8: Fig. S2D). We note that in the LUAD and LUSC datasets both PD-1 and CD80 expression correlated with IRF1, but not with IRF9, indicating the participation of IFN-γ signaling in both ligand and receptor expression during the development of immune resistance (Additional file 4: Table S3). RGMB expression was not correlated with IRF1 and IRF9 expression (Additional file 4: Table S3).
We next addressed how PD-L1 and PD-L2 expression correlated genes converged with three previously described gene signatures, IFN-γ, expanded immune, and T cell inflamed, used for predicting pembrolizumab responders (Additional file 9: Table S7) [40]. The genes belonging to these gene signatures were largely included among the PD-L2 expression correlated genes in the LUAD and LUSC datasets (Fig. 3a). In Lung_NSC the overlap was less evident as could be expected due to the lack of immune cells and cytokine stimulation in these in vitro grown cancer cells (Fig. 3a). The overlap between the three gene signatures and PD-L1 expression correlated genes was less pronounced than that for PD-L2 expression correlated genes (Fig. 3a). We also examined how previously described gene signatures representing immune cell types converged with PD-L1 and PD-L2 expression correlated genes in the LUAD, LUSC, and Lung_NSC datasets (Fig. 3a) (Additional file 9: Table S7) [21]. Again, the overlap between these gene signatures and PD-L1 expression correlated genes was less pronounced than that for PD-L2 expression correlated genes (Fig. 3a).
We analyzed the potential of PD-L1 and PD-L2 expression correlated genes to represent gene signatures with potential biomarker implications for immunotherapy with PD-1/PD-L1 axis blockade by searching for genes whose expression were correlated across the datasets Lung_NSC, LUAD, and LUSC. Notably, only PD-L1 and PD-L2 themselves remained after selecting for expression correlated genes for both PD-L1 and PD-L2 across all three datasets (Fig. 3b). PD-L2, plasminogen receptor with a c-terminal lysine (PLGRKT), and apolipoprotein L6 (APOL6) expression correlated with PD-L1 across the Lung_NSC, LUAD and LUSC datasets (Fig. 3c and Additional file 10: Table S8). For PD-L2, the three genes PD-L1, Transmembrane protein 106A (TMEM106A), and Tripartite motif containing 22 (TRIM22) expression correlated across the Lung_NSC, LUAD, and LUSC datasets (Fig. 3d and Additional file 10: Table S8). Notably, these genes also represented the only PD-L1 and PD-L2 expression correlated genes across the Lung_NSC and LUSC datasets (Fig. 3c, d and Additional file 10: Table S8). The expression correlation between the identified genes, PD-L1, and PD-L2 in individual LUAD and LUSC tumor samples is illustrated in Additional file 11: Fig. S3A and B. The examination of available LUSC and LUAD PD-L1 protein expression data also revealed expression correlation with PD-L1 protein except for APOL6 mRNA (Additional file 11: Fig. S3A and B). Next, we focused on PD-L1 expression correlated genes across the Lung_NSC and LUAD datasets (Additional file 10: Table S8). Here we identified 49 genes (Fig. 3c and Additional file 10: Table S8) and their expression correlation with PD-L1 for individual LUAD tumors is illustrated in Additional file 11: Fig. S3C. The mRNA expression for most of these 49 genes correlated with the expression of PD-L1 protein (Additional file 11: Fig. S3C, upper panel). For PD-L2, 26 genes expression correlated across the Lung_NSC and LUAD datasets (Fig. 3d and Additional file 10: Table S8). The expression correlation of these 26 genes with PD-L2 for individual LUAD tumor samples is illustrated in Additional file 11: Fig. S3C. The mRNA expression for most of these 26 genes correlated with the expression of PD-L1 protein (Additional file 11: Fig. S3C, lower panel).
In LUSC, PD-L1 expression correlated genes are enriched for being Chr9p24 localized
We noted fewer number of genes with PD-L1 expression correlation (n = 26) as compared to the number of genes with PD-L2 expression correlation (n = 326) in LUSC (Fig. 2e). GSEA for chromosomal localization of the 26 PD-L1 expression correlated genes in LUSC identified enrichment for localization to the Chr9p24 region for 13/26 genes (Fig. 4a, b and Additional file 12: Table S9). Chr9p24 enrichment was also present among the PD-L1 expression correlated genes in LUAD (6/257 genes), but not among the PD-L1 expression correlated genes in Lung_NSC (Fig. 4b and Additional file 12: Table S9). GSEA also identified enrichment for Chr9p24 localization among the PD-L2 expression correlated genes in LUSC, but for a lower fraction of genes (6/326) than observed for PD-L1 (13/26) (Fig. 4b and Additional file 12: Table S9). Enrichment for localization to additional genomic regions, e.g., Chr12p13 and Chr19q13, was also present among the PD-L1 and PD-L2 expression correlated genes in LUAD and LUSC (Additional file 12: Table S9). In LUSC, an expression correlation between PD-L1 protein and 7/12 examinable Chr9p24-localized genes with PD-L1 expression correlation was observed (note the Chr9p24-localized gene UHRF was excluded from the analysis) (Fig. 4c). In LUAD, the expression correlation between Chr9p24-localized genes with PD-L1 expression correlation and PD-L1 protein was less pronounced (Fig. 4d). We next performed an unsorted hierarchical cluster analysis of a merged gene signature composed of genes representing the core IFN signaling, different immune cell types, and the PD-L1 expression correlated genes in LUSC with Chr9p24 localization (the genes and the merged gene signature is shown in Additional file 9: Table S7). We found that in the LUSC samples the PD-L1 expression correlated genes with Chr9p24 localization were preferentially clustering together (Additional file 13: Fig. S4A). Furthermore, this cluster of Chr9p24 localized genes included PD-L1 and PD-L2. This point to Chr9p24 localization is the major determinant for the PD-L1 and PD-L2 expression profiles in LUSC (Additional file 13: Fig. S4A). The same analysis in LUAD samples showed that PD-L1 and PD-L2 did not cluster with the other Chr9p24-localized genes, demonstrating that in LUAD, the Chr9p24 localization was less important for determining the PD-L1 and PD-L2 expression profile (Additional file 13: Fig. S4B). Altogether, we conclude that the expression analyses of Chr9p24-localized genes with PD-L1 expression correlation may possess relevance when searching for gene expression biomarkers that can predict responsiveness to immunotherapy with PD-1/PD-L1 axis blockade in lung squamous cell carcinoma.