PD-1 and PD-L1 correlated gene expression profiles and their association with clinical outcomes of breast cancer

Background Immunotherapies that targeting programmed cell death 1 (PD-1) and programmed death-ligand 1 (PD-L1) have obtained prominent success in breast cancer (BC). However, not all the patients benefit from the antibody therapy. This study aimed to identify PD-1/PD-L1 correlated genes and pathways as well as investigate their potential as prognostic marker in BC. Materials and methods By analysing transcriptional data of BC from TCGA, we identified PD-1 and PD-L1 correlated genes by WGCNA analysis and explored the biological process as well as pathways they enriched. Co-expression analysis were performed for PD-1/PD-L1 with immune infiltration and checkpoints. The prognostic value of PD-1 and PD-L1 were also investigated. Results PD-1 and PD-L1 expression showed significant difference in different molecular subtypes and stages. PD-1 correlated genes enriched in T cell activation, lymphocyte activation, leukocyte migration while PD-L1 correlated genes demonstrated enrichment including T cell apoptotic process, tolerance induction and cytolysis. Immune infiltration analysis suggested that PD-1 and PD-L1 were related with Neutrophils (r = 0.65, r = 0.48) and Fibroblasts (r = 0.59, r = 0.47). For immune checkpoints analysis, PD-1 was associated with HLA-A (r = 0.804) and INPP5D (r = 0.782) while PD-L1 correlated with CTLA4 (r = 0.843) and CD27 (r = 0.823). PD-1 was associated favorable survival of BC (HR = 0.67, P = 0.012) while PD-L1 did not demonstrate significant association with BC prognosis (HR = 0.85, P = 0.313). Conclusion PD-1 and PD-L1 correlated genes participated in biological process including T cell activation, lymphocyte activation, leukocyte migration, T cell apoptotic process, tolerance induction and cytolysis. PD-1/PD-L1 expression also demonstrated relation with immune infiltration and immune checkpoints. High PD-1 expression predicted better survival of breast cancer patients.


Background
As one of the most frequently occurred malignant tumors, breast cancer (BC) remains the leading cause of cancer-related death for females in many countries [1]. Breast cancer arises from multiple genetic factors, environmental alternations and their complicated interactions [2]. According to the status of biomarkers including estrogen receptor (ER), progesterone receptor (PR) and epidermal growth factor receptor 2 (HER2), patients with breast cancer were classified into groups of luminal A, luminal B, HER-2 positive and triple negative [3]. It has been accepted that different groups of BC patients benefit from corresponding treatment strategy of chemical and hormonal therapy [4].
Although certain therapeutic combinations have been used as standard treatment in clinical management of BC, some BC patients still could not get satisfactory clinical outcomes [5]. The different outcomes of BC patients indicated that other critical factors also determine the final therapeutic effect such as the immune status of the cells [6]. It is well-accepted that immune escape of tumor cells and aberrant human immune surveillance play essential role in carcinogenesis, progression and metastasis of various types of cancer [7]. As for immune escape of cancer cells, the identification of PD-1 (programmed death 1) and PD-L1 (programmed death-ligand 1) axis was one of the most encouraging finding of cancer therapy in recent years [8]. Serving as an immune checkpoint in tumor microenvironment, the antibodies of PD-1/ PD-L1 has shown prominent effect in a large number of cancer types [9].
Previously, PD-L1 expression has been reported to be associated with worse prognosis of triple negative breast cancer patients, which counteract effect of tumorinfiltrating lymphocytes (TILs) [10]. Another study of HER2 + invasive BC patients indicated a positive correlation of PD-L1 expression and CD8+ T cells with favorable clinical outcomes [11]. One research of 1318 BC patients in European suggested that PD-1 positive immune cells in triple negative breast cancer correlated with longer disease-free survival, and tumor-infiltrating lymphocytes (TILs) density was remarkably related with PD-1 and PD-L1 expression in immune cells [12]. Most of the findings implied that medicine concerning PD-1 and PD-L1 immune checkpoint might become novel therapeutic strategies for breast cancer.
Although the critical role of immune checkpoint PD-1 and PD-L1 have been widely reported in a number of malignant tumors, the underlying regulating mechanisms in breast cancer is still unclear. In this study, we performed comprehensive analysis of gene expression profiles related to PD-1 and PD-L1 in breast cancer using transcriptome data from TCGA. The correlation of PD-1 and PD-L1 with other immune biomarkers and immune cells infiltration were revealed. Furthermore, the effect of PD-1 and PD-L1 on clinical outcomes of BC was explored to determine their potential as biomarkers for BC patients prognosis.

Analyzed datasets
The RNA sequencing and clinical data of breast cancer patients in TCGA datasets were downloaded from UCSC XENA (https ://xena.ucsc.edu/). The level of gene expression was measured as Transcripts per million reads (TPM). Clinical data included the histological type, molecular Type, cancer stage, recurrence event and survival information. The relationship between PD-1/PD-L1 expression and the clinical data were investigated.

Co-expression gene and enrichment analysis
Weighted correlation network analysis (WGCNA) is an algorithm for finding genetic interactions in a weighted manner. Co-expressed genes obtained by WGCNA analysis will be more accurate. Using WGCNA analysis, we searched the co-expressed genes for PD-1 and PD-L1. As genes with little variation in expression usually represent noise, the most variant genes were filtered for network construction. Gene variabilities were measured by median absolute deviation (MAD). If the interaction gene was more than 200, we used the interaction degree to search the top 200 gene as the interaction genes. Clusterprofiler is a R package for enrichment analysis. Using clusterprofier, we used biological process in the Gene Ontology (GO) to analyze the interacted genes [13]. Since the results of the enrichment analysis contain many similar results, we further concentrated the results of the enrichment analysis.

Relationship between immune factors and PD-1/PD-L1
A variety of studies have confirmed that immune infiltration and all aspects of the tumor are related. MCP-counter is available R package to estimate the sample immune infiltration. From a gene expression matrix, it produces for each sample an abundance score for CD3+ T cells, CD8+ T cells, cytotoxic lymphocytes, NK cells, B lymphocytes, cells originating from monocytes (monocytic lineage), myeloid dendritic cells, neutrophils, as well as endothelial cells and fibroblasts. Then we used correlation analysis to evaluate the correlation between PD-1/ PD-L1 and immune infiltration. Immunological checkpoints serve as the primary site for detecting immune status, and we also evaluated the relationship between PD-1 and PD-L1 and immune checkpoints.

Statistical analysis
In this study, statistical analysis was mainly performed by using R language (https ://www.r-proje ct.org/) with several publicly available packages. Rank sum test was used to evaluate the expression difference of PD-1/PD-L1 in different groups. Spearman correlation analysis was used to explore the correlation between PD-1/PD-L1 expression and immune infiltration and immune checkpoints. Survival curve was generated by Kaplan-Meier method based on log-rank test. Other Figures were generated by several R packages, such as pheatmap, circlize, and corrplot. All multiple tests were corrected by the BH method. A probability value P < 0.05 was considered to be significant in this study.

PD-1/PD-L1 expression status in different clinical subgroups
Using TCGA datasets, we analyzed the PD-1/PD-L1 expression in different groups according to the clinical data. As shown in Fig. 1a, both PD-1 and PD-L1 expression showed significant difference in the three molecular subtypes (P < 0.001 and P = 0.047, respectively). Luminal and Basal-like subtype show significant difference in PD-1 and PD-L1 expression. Moreover, the expression of PD-L1 differs among different stages (P = 0.032), while PD-1 did not show any difference (P = 0.536). In addition, both PD-1 and PD-L1 expression correlated with the recurrence event of BC patients (P = 0.017, P = 0.015, respectively). We also analyzed the relation of PD-1/ PD-L1 expression with clinical data including therapy, histological subtype, ER, PR and HER-2 status (Table 1). Finally, both PD-1 and PD-L1 were correlated with ER, PR and clinical therapy, indicating the probable implication of PD-1/PD-L1 in clinical outcome.

Co-expression analysis of genes associated with PD-1 and PD-L1
Using WGCNA, we analyzed the co-expression gene associated with PD-1 and PD-L1. The connectivity among a b Fig. 1 Differences in PD-1 and PD-L1 expression between different clinicopathological information. a PD-1 and PD-L1 expression in molecular type. b PD-1 and PD-L1 expression in stage genes was a scale-free network distribution if the value of soft thresholding power β equals to 3 (Fig. 2a). Altogether 21 module was obtained according to WGCNA analysis (Fig. 2b). Among these modules, PD-1 belonged to pink module while PD-L1 belonged to thistle 1 module. We finally got 1065 genes that interacted with PD-1 and 99 PD-L1 correlated genes. Then we selected the top 200 gene associated with PD-1 and all of the 99 PD-L1 related   Table 2).

PD-1/PD-L1 expression and immune infiltration
Using Microenvironment Cell Populations-counter, we evaluated the profiles of immune infiltration among various subtypes and stages breast cancer (Fig. 3a). Additionally, the associations of PD-1 and PD-L1 with immune cell populations according to the transcriptomic data were analyzed. The results indicated that PD-1 and PD-L1 were mainly related with Neutrophils (r = 0.65, r = 0.48) and Fibroblasts (r = 0.59, r = 0.47) (Fig. 3b).

Discussion
It is well-accepted that tumor and its microenvironment is a complex unit to complete the aggressive growth and metastasis of cancer cells [14][15][16]. Although a number of chemical and radical therapy have been used in clinical practice of cancer treatment, the immune system of human are still believed to be the most fundamental and effective weapon against cancer [17]. Programmed cell death 1 (PD-1) and corresponding ligand (PD-L1) are one of the most critical biological suppressors of cytotoxic immune reaction, the antibody of which has been authorized by FDA in a unprecedented fast term [18,19]. However, the inhibitor of PD-1/PD-L1 were not effective in every individual, which require the comprehensive understanding of specific mechanisms underlying PD-1/ PD-L1 regulation in carcinogenesis [20].
Using data from TCGA, we first unraveled the expression status of PD-1/PD-L1 in different subtypes and clinical stages of BC patients. PD-L1 expression was decreased in stage III-IV compared with stage I-II. On the basis of molecular classification, basal-like BC subtype showed highest expression of PD-1 and ERBB2 subtype BC had highest PD-L1 expression, which suggest that different subtypes possess various PD-1/PD-L1 status. In addition, PD-1 and PD-L1 expression correlated with the recurrence of BC patients, which also confirm the critical role of PD-1/PD-L1 immune checkpoint in BC progression.
We subsequently identified co-expressed genes of PD-1 and PD-L1 by means of WGCNA. A total of 1065 genes correlated with PD-1 and 99 PD-L1 correlated genes were screened. The GO enrichment analysis of PD-1 correlated genes suggested biological process of T cell activation, regulation of lymphocyte activation, regulation of T cell activation and leukocyte migration.  the tumor tissue, leading to an adaptive upregulation of PD-L1 and tumor protection against immune destruction [21]. PD-L1 is expressed by antigen-presenting cells and results in T cell inactivation by interaction with PD-1 on T-cells. It is therefore of great importance to clarify the biological process and pathways we identified for PD-1/ PDL1 related genes, by which immunotherapy with PD-1-and PD-L1-targeted monoclonal antibodies might dramatically change the therapeutic and prognostic landscape for cancer. Immune infiltration of breast cancer determine the immune activation of tumor microenvironment and is related with clinical outcome of patients. According to the correlation analysis of PD-1 and PD-L1 with immune cell populations, PD-1 and PD-L1 were mainly related with Neutrophils and Fibroblasts. It has been reported that PD-1 protein expression significantly correlated with higher TIL abundance, Ki-67 index, basallike subtypes, and distant metastasis of triple-negative breast cancer (TNBC) [22]. The interaction of PD-1/ PD-L1 immune checkpoint with immune infiltration is a promising research direction in the future. The immune checkpoints of CD28, CD80, CD86, CTLA4, INPP5D, INPPL1, CD58, CD27, CD70, HLA-A, CD74 were also analyzed in relation to PD-1/PD-L1 expression. Finally, PD-1 was found to be mainly associated with HLA-A and INPP5D while PD-L1 correlated with CTLA4 and CD27.
Cytotoxic T-lymphocyte-associated protein 4 (CTLA4) belongs to immunoglobulin superfamily and encodes a protein transmitting inhibitory signal to T cells, the antibody of which also demonstrated favorable outcome in clinical management [23]. The combination usage of PD-1/PD-L1 with these immune checkpoints might become novel therapeutic targets in the future [24].
Survival analysis of PD-1/PD-L1 expression status demonstrated that the expression of PD-1 was associated favorable survival of breast cancer patients while PD-L1 did not suggest significant association with BC prognosis. In a study of 195 triple-negative breast cancer individuals, PD-1 was found to be significantly related with better disease free survival and overall survival. The results of this study also suggested that PD-1 protein expression in TILs, but not PD-L1 in tumor cells, predicted better prognosis in TNBC [22]. Researchers demonstrated that high expression of PDL1 are associated with favorable clinical outcome in 127 primary breast cancer [25]. Another study of HER2+ invasive BC patients indicated a positive correlation of PD-L1 expression and CD8+ T cells with favorable clinical outcomes [11]. On contrary, PD-L1 expression has been reported to be associated with worse prognosis of triple negative breast cancer patients, which counteract effect of tumorinfiltrating lymphocytes [10]. Response of BC patients to neoadjuvant therapy and survival outcome indicated that PD-L1 predicted better rate of pathological complete response (pCR) [26]. In addition, PD-1/PD-L1 correlated genes such as CD5, CD74, CD96 and CD226 were also related with prognosis of BC patients. The combination of these immune markers with PD-1/PD-L1 might improve the prediction and management of BC patients in the future. Until now, whether PD-1/PL-L1 predict prognosis in BC patients was still in debate. The specific association of PD-1/PD-L1 with BC prognosis still required further large-scale investigations to confirm, which would benefit the potential of PD-1/PD-L1 as prognostic biomarker.

Conclusion
In this study, we reported the expression status of PD-1 and PD-L1 in different subtypes of breast cancer. The PD-1/PDL1 correlated gene profiles were described, the enrichment analysis of which focus on biological process including T cell activation, lymphocyte activation, leukocyte migration, T cell apoptotic process, tolerance induction and cytolysis. PD-1/PD-L1 expression also demonstrated relation with immune cells infiltration and multiple immune checkpoints. High PD-1 expression was associated with better survival of breast cancer patients.