PD-L1 and PD-L2 expression correlated genes in non-small-cell lung cancer

Background Programmed cell death ligand-1 (PD-L1) and ligand-2 (PD-L2) interaction with programmed cell death protein-1 (PD-1) represent an immune-inhibiting checkpoint mediating immune evasion and is, accordingly, an important target for blockade-based immunotherapy in cancer. In non-small-cell lung cancer (NSCLC), improved understanding of PD-1 checkpoint blockade-responsive biology and identification of biomarkers for prediction of a clinical response to immunotherapy is warranted. Thus, in the present study, we systematically described PD-L1 and PD-L2 expression correlated genes in NSCLC. Methods We performed comparative retrospective analyses to identify PD-L1 and PD-L2 mRNA expression correlated genes in NSCLC. For this, we examined available datasets from the cancer cell line encyclopedia (CCLE) project lung non-small-cell (Lung_NSC) and the cancer genome atlas (TCGA) projects lung adenocarcinoma (LUAD) and squamous cell carcinoma (LUSC). Results Analysis of the CCLE dataset Lung_NSC identified expression correlation between PD-L1 and PD-L2. Moreover, we identified expression correlation between 489 genes and PD-L1, 191 genes and PD-L2, and 111 genes for both. PD-L1 and PD-L2 also expression correlated in TCGA datasets LUAD and LUSC. In LUAD, we identified expression correlation between 257 genes and PD-L1, 914 genes and PD-L2, and 211 genes for both. In LUSC, we identified expression correlation between 26 genes and PD-L1, 326 genes and PD-L2, and 13 genes for both. Only a few genes expression correlated with PD-L1 and PD-L2 across the CCLE and TCGA datasets. Expression of Interferon signaling-involved genes converged in particular with the expression correlated genes for PD-L1 in Lung_NSC, for PD-L2 in LUSC, and for both PD-L1 and PD-L2 in LUAD. In LUSC, PD-L1, and to a lesser extent PD-L2, expression correlated with chromosome 9p24 localized genes, indicating a chromosome 9p24 topologically associated domain as an important driver of in particular LUSC PD-L1 expression. Expression correlation analyses of the PD-L1 and PD-L2 receptors programmed cell death protein-1 (PD-1), Cluster of differentiation 80 (CD80), and Repulsive guidance molecule B (RGMB) showed that PD-1 and CD80 expression correlated with both PD-L1 and PD-L2 in LUAD. CD80 expression correlated with PD-L2 in LUSC. Conclusions We present gene signatures associated with PD-L1 and PD-L2 mRNA expression in NSCLC which could possess importance in relation to understand PD-1 checkpoint blockade-responsive biology and development of gene signature based biomarkers for predicting clinical responses to immunotherapy. Electronic supplementary material The online version of this article (10.1186/s40880-019-0376-6) contains supplementary material, which is available to authorized users.

Tumor PD-L1 expression in NSCLC has been shown to be associated with poor prognosis in two meta-analyses [23,24]. Moreover, meta-analyses of clinical trials have illustrated an overall positive response to immunotherapy with PD-L1/PD-1 axis blockade in NSCLC patients and positive correlation between the PD-L1 tumor proportion score (TPS) (the percentage of PD-L1 expressing tumor cells determined by immunohistochemistry) and the response rate [25,26]. In the CheckMate 017 and CheckMate 057 trials corresponding to squamous and non-squamous NSCLC patients, respectively, patients receiving the PD-1 antibody, nivolumab, as a second-line of treatment showed a response rate of 20% regardless of PD-L1 TPS [4,5]. Analysis of subgroups defined by PD-L1 expression showed that PD-L1 TPS ≥ 1%, ≥ 5%, and ≥ 10% were associated with improved response (See figure on next page.) Fig. 1 PD-L1 and PD-L2 expression correlation analyses in 114 NSCLC cell lines from the CCLE dataset (Lung_NSC). a Overview of key factors interacting with PD-L1 and PD-L2 in NSCLC. PD-L1 and PD-L2 are ligands of PD-1 and the interaction of PD-L1 or PD-L2 with PD-1 results in T cell suppression. Moreover, PD-L1 can interact with CD80 on activated T cells and inhibit T cell activity. PD-L2 has a second receptor, RGMB, and this interaction can activate T cells. b Expression correlation analysis of PD-L1 and PD-L2 in the CCLE dataset (Lung_NSC, n = 114) using Broad portal analysis of RNA sequencing-based gene expression data. Spearman and Pearson correlation coefficients r and the corresponding P value are shown. c Gene expression correlation between PD-1, CD80, RGMB, PD-L1, and PD-L2 in CCLE dataset (Lung_NSC, n = 114) Spearman and Pearson correlation coefficients r and the corresponding P values are shown. The horizontal lines show the pairwise comparison of genes and line thickness is proportional to the degree of expression correlation. Dashed lines indicate negative expression correlation. Significant expression correlation is indicated in red. d Number of genes expression correlated with PD-L1, PD-L2, and both in CCLE dataset (Lung_NSC, n = 114). For all panels, the criteria for significant expression correlation were Pearson correlation coefficient r ≥ 0.3 or ≤ − 0.3, Spearman correlation coefficient r ≥ 0.4 or ≤ − 0.4, and P values < 0.05. The analysis was performed using the GenomicScape portal. BMP bone morphogenic protein, BMPR I BMP receptor type 1, BMPR II BMP receptor type 2, C2-type constant 2-type, CD80 cluster of differentiation 80, IFN interferon, Ig-like Immunoglobulin-like, PD-1 programmed cell death protein-1, PD-L1 programmed cell death ligand-1, PD-L2 programmed cell death ligand-2, RGMB repulsive guidance molecule B, SHP-1/2 Src homology 2 domain-containing protein tyrosine phosphatase 1 and 2, Pe Pearson, r correlation coefficient, Sp Spearman, V-type variable-type rates (31%, 36%, and 37%, respectively) in the Check-Mate 057, whereas, an association between PD-L1 TPS and response rate was not evident from the CheckMate 017 [4,5]. In the KEYNOTE-001 and KEYNOTE-24 trials, which were the first trials to use first-and secondline therapies with PD-1 antibody, pembrolizumab, in advanced NSCLC patients, PD-L1 TPS ≥ 50% resulted in response rates of 52% and 45%, respectively [27][28][29]. The PD-L1 antibodies, atezolizumab, durvalumab, and avelumab, also showed promising response rates as immunotherapy in NSCLC even with low PD-L1 TPS [6,30]. Unlike the PD-1 antibodies, the PD-L1 antibodies also block the PD-L1/CD80 axis but leave the PD-L2/PD-1 axis unaffected. The response rate in NSCLC was similar for PD-1 and PD-L1 antibody immunotherapy but less immune-related adverse effects were observed from the latter [31]. The efficient stratification of immunotherapy responders and non-responders among NSCLC patients remains a major challenge [2,3,6,32]. Given the potential overlapping involvement of PD-L1 and PD-L2 for defining immune checkpoint inhibition, we performed this retrospective expression correlation study to identify genes in NSCLC having expression profiles closely associated with the PD-L1 and PD-L2 expression profiles. The presented study illustrates novel aspects of PD-L1 and PD-L2 regulation, with potential biological relevance, as well as relevance for immunotherapy response stratification.

Datasets
The following datasets were used for retrospective gene expression analyses. The cancer cell line encyclopedia (CCLE) project dataset is a compilation of gene expression data from human cancer cell lines [33]. The CCLE subset, lung_non-small-cell (Lung_NSC), contains gene expression data from NSCLC cell lines, were obtained from the CCLE portal (https ://porta ls.broad insti tute.org/ ccle/home). The Cancer Genome Atlas (TCGA) project dataset is a comprehensive atlas of gene expression and gene regulation across human cancers [34]. The TCGA data subsets, lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC), contain gene expression data from dissected NSCLC tumors, were obtained from the TCGA data portal (https ://tcga-data.nci.nih. gov). The Genotype-Tissue Expression (GTEx) project dataset is a comprehensive atlas database of gene expression in human tissues [35]. The GTEx lung subset, contains RNA sequencing-based expression data from normal lung tissue, were obtained from the GTEx portal (http://www.gtexp ortal .org/home/).

Broad portal analysis
The Broad portal (https ://porta ls.broad insti tute.org/ccle) was used to generate scatter plots from CCLE Lung_NSC expression data for PD-L1 and PD-L2 using Illumina RNA sequencing data (processed with normalization followed by Log2 transformation).

cBioPortal analysis
TCGA RNA sequencing data (in format V2 RSEM Z-score) were in the cBioPortal for Cancer Genomics interface (http://www.cbiop ortal .org/) used to generate gene lists composed of expression correlated genes for PD-L1, PD-L2, IRF1, and IRF9. Correlation coefficients were calculated before the Log2 transformation of RNA sequencing data. Log2 transformed TCGA RNA sequencing data were used for heat map analyses and hierarchical clustering. Finally, the plot function in the cBioPortal interface was used to calculate correlation coefficients between PD-L1 protein expression (from reverse-phase protein arrays) and mRNA expression of identified PD-L1 and PD-L2 expression correlated genes.

Xena portal analysis
UCSC Xena browser (https ://xenab rowse r.net/) was used to calculate pairwise expression correlation using the GTEx lung and TCGA RNA sequencing data (in format Log2[normalized_counts+1]). Expression correlation coefficients and corresponding P values were calculated before Log2 transformation using the browser http:// gepia .cance r-pku.cn/ which analyses data present in the Xena portal.

Gene set enrichment analysis (GSEA)
The molecular signatures database (MSigDB) is a collection of annotated gene sets [36,37]. GSEA was used to compute overlaps between gene sets representing expression correlated genes for PD-L1, PD-L2, Interferon response factor 1 (IRF1), and Interferon response factor 9 (IRF9) and available gene sets in MSigDB (http://softw are.broad insti tute.org/gsea/). We analyzed the MSigDB gene set collection named hallmark gene sets, which summarized and represented specific well-defined biological states or processes [36,37]. Moreover, we analyzed the MSigDB gene set collection named positional gene sets, corresponding to each human chromosome and each cytogenetic band that contained at least one gene. Values in the tables presenting GSEA represented the number of genes in the respective MSigDB hallmark gene sets (K), the number of genes in the examined PD-L1, PD-L2, IRF1, and IRF9 expression correlation gene signatures overlapping with the MSigDB hallmark gene set (k), P value representing the hypergeometric distribution of overlapping genes, and the false discovery rate (FDR) q value after correction for multiple hypothesis testing. Enriched gene signatures with q value < 1E−03 were considered significant. The list of the specific gene sets analyzed and their sources are available in supplementary tables.

Statistical analysis
To decrease data overfitting, the criteria used for significant expression correlation between two genes was a Pearson correlation coefficient r ≥ 0.3 or ≤ − 0.3, a Spearman correlation coefficient r ≥ 0.4 or ≤ − 0.4, and all corresponding P values < 0.05. Notably, cBioPortal and Xena portal analyses included gene expression data for different numbers of TCGA LUAD and LUSC patients resulting in minor differences in obtained correlation coefficients and resulting statistics.

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-γ  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).
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.

Discussion
Comparison of PD-L1 and PD-L2 gene expression results from the CCLE Lung_NSC cell line dataset and the TCGA LUAD and LUSC tumor datasets include some Expression analyses using the CCLE dataset identified gene expression in various cancer cell genetic backgrounds each with a unique and relative homogenous profile in terms of cancer driver and passenger mutations. However, the actual cell propagation conditions only partially mimic the context of a tumor in vivo. Gene expression analyses using the TCGA datasets (LUAD and LUSC) corresponded to incisional and core tumor biopsies, which mimic tumor cells more naturally due to stimulation of the tumor cells with factors, such as cytokines, derived from the complex tumor microenvironment. However, these biopsies may not be representative for the entire tumor burden due to tumor heterogeneity and they may as well display various degrees of tumor-infiltrating immune cells and connective tissue. With such dataset differences in mind, we still found it surprising that only PD-L2, APOL6, and PLGRKT expression were correlated with PD-L1 across the CCLE Lung_NSC and TCGA (LUAD and LUSC) datasets. Similar, only the three genes PD-L1, TRIM22, and TMEM106A expression correlated with PD-L2 across the three datasets. As a result, PD-L2 expression should in principle be among the best indirect predictors to deduce PD-L1 mRNA expression level, and vice versa. This supports the use of PD-L2 expression analyses, beyond PD-L1 expression analyses, to predict responsiveness to immunotherapy with PD-1/PD-L1 axis blockade in NSCLC.
Examining PD-L1 and PD-L2 gene expression in the NSCLC datasets separately revealed other aspects of expression regulation. In the CCLE dataset Lung_NSC, the majority of the PD-L2 expression correlated genes were also expression correlated with PD-L1, 111/191 genes, but not vice versa, 111/489 genes. In the TCGA dataset LUAD, the opposite was observed, 211/914 genes and 211/257 genes, respectively. In the TCGA dataset LUSC, we observed a low number of PD-L1 expression correlated genes (n = 26) relative to the number of genes expression correlated with PD-L2 (n = 326). Half of these PD-L1 expression correlated genes, 13/26, expression correlated with PD-L2. Thus, gene signatures representing PD-L1 and PD-L2 expression correlated genes seem to delineate common, as well as distinct, aspects of PD-L1 and PD-L2 regulation; depending on the actual NSCLC datasets analyzed.
With previous studies showing that both PD-L1 and PD-L2 genes are regulated by IFN signaling, we focused on IFN signaling as an important driving mechanism for the expression of PD-L1 and PD-L2 [21,40]. In the CCLE dataset (Lung_NSC) the expression correlation between IFN signaling signature genes and PD-L1 and PD-L2 differed strikingly to the expression correlation observed in the TCGA datasets (LUAD and LUSC). IRF9 expression correlated with PD-L1 in Lung_NSC, but not in the LUAD and LUSC datasets. However, IRF1 expression correlated with PD-L1 and PD-L2 expression in the LUAD and LUSC datasets indicated a shift from IFN-α involved PD-L1 expression regulation in CCLE (Lung_ NSC) cell lines toward IFN-γ involved PD-L1 and PD-L2 expression regulation in TCGA (LUAD and LUSC) tumors. IRF1 expression correlation with PD-L1 and PD-L2 also goes beyond lung cancer [21]. Furthermore, our presented expression correlation data indicated a stronger correlation of PD-L2 expression, relative to PD-L1 expression, with active IFN signaling in LUAD and LUSC tumors.
We found it intriguing that PD-L1 expression correlated with only a few other genes in the TCGA LUSC dataset, and that a prominent number of these genes were localized to Chr9p24. Of the genes located to Chr9p24 and expression correlated with PD-L1, only PD-L2 and JAK2 were associated with IFN-γ signaling. JAK2 is essential for intracellular signaling upon exposure of cancer cells to IFN-γ produced by T cells upon antigen recognition. In cancer, loss of function mutations in JAK2 are relatively often observed and can in principle result in decreased antigen presentation, lack of T cell infiltrates, tumor cell resistance to the anti-proliferative effects of IFN-γ, and can prevent adaptive PD-L1 and PD-L2 expression upon IFN-γ exposure [41,42]. One consequence could be intrinsic and/or acquired resistance to immunotherapy with PD-1/PD-L1 axis blockade [41,42]. In LUSC, Chr9p21.3 allelic deletions of the tumor suppressor gene Cyclin-dependent kinase inhibitor 2A (CDKN2A) can, in addition, include co-deletion of JAK2, PD-L1, PD-L2 [43]. Moreover, co-amplification of JAK2 and PD-L1 is observed in NSCLC and patients with JAK2 and PD-L1 co-amplification has a poorer overall survival and recurrence-free survival compared to patients with normal PD-L1 copy number [44]. Copy number variation (CNV) is a characteristic of NSCLC and drives the expression profile of underlying genes [45][46][47][48][49]. However, we found no evidence that Chr9p CNVs should be the driving mechanism for the distinct PD-L1 versus PD-L2 expression correlation profile in LUSC, as such CNVs should affect PD-L1 and PD-L2 expression simultaneous. The architecture of the human genome is a major determinant for controlling gene expression [50,51]. First, local enhancer-promoter interaction at the single gene level. Second, the genome is organized into functionally distinct topologically associated domains (TADs), from 100 kilo-bases to 5 mega-bases, encompassing multiple genes and regulatory elements [50,52,53]. TAD boundaries act as insulators by preventing communication between regulatory elements on either side, and TADs are mostly invariant between different cell types, while loops within TADs (sub-TADs) are more tissue-specific [50,52,53]. Finally, interactions between TADs organize the genome into mega-base compartments comprising open gene-rich chromatin or closed gene-poor chromatin [54]. In LUSC, the expression correlation between PD-L1 and Chr9p24 genes indicate the existence of a specific TAD architecture making PD-L1, but not PD-L2, expression largely dependent on enhancers shared with other Chr9p24 localized genes within this sub-TAD. One consequence is the possibility of PD-L1 expression in a non-IFN stimulated tumor environment or despite genetic defects in the IFN signaling cascade. Thus, gene expression profiling of Chr9p24 could have potential as an IFN-γ signaling independent predictive biomarker for PD-L1 expression with relevance for predicting responders to immunotherapy with PD-1/PD-L1 axis blockade in lung squamous cell carcinoma.
The limitations of this study are as follows. This was a retrospective study focusing on available gene expression data from TCGA NSCLC tumors and CCLE NSCLC cell lines. Gene expression data originated from different experimental platforms and raw data initially processed with different normalization strategies. Consequences related to potential data overfitting and lack of optimal data harmonization could negatively impair the validity of the subsequent comparisons. Nevertheless, the identified gene-signatures could be an important basis from which individual genes can be selected for further validation at mRNA and protein level in subsequent expression based analyses. In this line, the inclusion of further clinical data will be valuable to substantiate the usability of expression analyses of the hereby identified genes and gene-signatures as potential new clinical relevant response and resistance biomarkers for immunotherapy with PD-1/PD-L1 axis blockade in NSCLC.