Butyrophilin-like 9 expression is associated with outcome in lung adenocarcinoma
BMC Cancer volume 21, Article number: 1096 (2021)
Lung adenocarcinoma (LUAD) is the most prevalent non-small cell lung cancer (NSCLC). Patients with LUAD have a poor 5-year survival rate. The use of immune checkpoint inhibitors (ICIs) for the treatment of LUAD has been on the rise in the past decade. This study explored the prognostic role of butyrophilin-like 9 (BTNL9) in LUAD.
Gene expression profile of buytrophilins (BTNs) was determined using the GEPIA database. The effect of BTNL9 on the survival of LUAD patients was assessed using Kaplan-Meier plotter and OncoLnc. Correlation between BTNL9 expression and tumor-infiltrating immune cells (TILs) was explored using TIMER and GEPIA databases. Further, the relationship between BTNL9 expression and drug response was evaluated using CARE. Besides, construction and evaluation of nomogram based on BTNL9 expression and TNM stage.
BTNL9 expression was downregulated in LUAD and was associated with a poor probability of 1, 3, 5-years overall survival (OS). In addition, BTNL9 expression was regulated at epigenetic and post-transcriptional modification levels. Moreover, BTNL9 expression was significantly positively correlated with ImmuneScore and ESTIMATEScore. Furthermore, BTNL9 expression was positively associated with infiltration levels of B cells, CD4+ T cells, and macrophages. Kaplan-Meier analysis showed that BTNL9 expression in B cells and dendritic cells (DCs) was significantly associated with OS. BTNL9 expression was significantly positively correlated with CARE scores.
These findings show that BTNL9 is a potential prognostic biomarker for LUAD. Low BTNL9 expression levels associated with low infiltration levels of naïve B cells, and DCs in the tumor microenvironment are unfavorable for OS in LUAD patients.
Lung cancer is the most common cancer and the leading cause of cancer-related deaths globally and in China [1, 2]. Although the 5-year survival rate has increased over the past four decades, the OS is poor (5.6–20.6%) . Immunotherapies have significantly improved cancer treatment during the past decade. For example, pembrolizumab used to treat naive advanced non–small-cell lung cancer (NSCLC) shows a 5-year survival rate of 23.2 and 29.6% in patients with a PD-L1 tumor proportion score ≥ of 50% . Immune checkpoint inhibitors (ICIs) block immune checkpoint signaling, thus alleviating antitumor immunity, and significantly improving five year-OS in NSCLC.
PD-1/PD-L1 is the most widely used ICIs, whereas other immune checkpoints, such as LAG-3, TIGIT, TIM-3, and CTLA-4, are currently under development . However, the biological role of immune checkpoint buytrophilins (BTNs)  in the regulation of NSCLC physiology and its underlying molecular mechanism remains to be fully elucidated. BTNs, including butyrophilin (BTN) and butyrophilin-like (BTLN), are related to the B7 family of co-stimulatory molecules. This family plays a significant role in T cell suppression, regulating epithelial cell and T cell interplays . Human BTN genes are located in the MHC class I domain of the short arm of chromosome 6 (6p22.1). Human BTN genes are grouped into three subfamilies, which form phylogenetically related groups, including BTN1, BTN2, and BTN3. BTN1A1 belongs to the BTN1 subfamily, BTN2A1, BTN2A2, and BTN2A3 (BTN2A3P) belong to the BTN2 subfamily, whereas BTN3A1, BTN3A2, and BTN3A3 belong to the BTN3 subfamily. Moreover, butyrophilin-like proteins (BTNL: BTNL2, BTNL3, BTNL8, BTNL9, and BTNL10) and SKINT-like (SKINTL) are classified in the family of BTNs [7,8,9].
Lung adenocarcinoma (LUAD) has been the most prevalent histopathological subtype of NSCLC in China since 2014 . In this study, we explored the relationship between the expression level of BTNs and LUAD prognosis. Significant survival-related BTNs were screened using GEPIA . Datasets used for analysis in this study were retrieved from Gene Expression Omnibus , TIMER , KM plotter , UALCAN , OncoLnc , Oncomine , TissGDB  databases. The findings from this study provide information on the relationship between immune checkpoint BTNL9 and tumor immune response. These findings show that BTNL9 can be used for the prognosis and development of immunotherapy for LUAD. A flow chart of the study design is shown in Fig. 1.
Determination of expression profiles of BTNs genes
Expression profiles of BTNs in LUAD were determined using the GEPIA database . Default settings were used with |Log2FC| > 1 and p-value Cutoff < 0.01 used as the cutoff criteria to determine differentially expressed genes. Log2(TPM + 1) was transformed for gene expression profile; jitter size was set at 0.4 for the plot. Notably, TCGA and GTEx datasets were included in the analysis. mRNA expression profile of BTNL9 in different tumor types was determined using the Oncomine database  using default settings with a P-value of 0.01, a fold change of 1.5, and a top 10% gene ranking.
Median gene expression was used as the cutoff point for survival analysis. Survival analysis of BTNL9 was performed using three online databases, including KM plotter , UALCAN , and OncoLnc databases . TissGDB is a tissue-specific gene annotation database in cancer . Forest plot were generated using Cox proportional hazard ratio (HR), and overall survival (OS) and relapse-free survival (RFS) of 28 cancer types were performed using 95% CI.
Estimation of infiltration level of immune cell type and correlation with BTNL9
TIMER is a comprehensive resource for systematically analyzing immune infiltration in diverse cancer types based on the TCGA dataset [13, 20]. The Gene_DE module from TIMER was used to calculate BTNL9 mRNA expression in cross-carcinoma (*: P-value < 0.05; **: P-value < 0.01; ***: P-value < 0.001). The expression profile of BTNL9 in LUAD and its correlation with six immune infiltration cells, including B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and DCs, were analyzed using Gene and Survival module in TIMER. Gene_Corr module was used to determine the correlation between BTNL9 expression and B and DC cells . The immune score and stromal score of each TCGA tumor sample were estimated using Sangerbox (http://sangerbox.com/Index).
Predicting binding of miRNA and lncRNA to BTNL9
miRMap , TargetScan , and miRWalk  were used to predict miRNAs that can bind to BTNL9. Predicted miRNAs obtained from the three databases were further verified using the starBase database .
Gene set enrichment analysis (GSEA) of BTNL9 high and low expression groups
Sangerbox is a tool developed by Hangzhou Mugu Technology Co., Ltd. GSEA was used to perform KEGG and HALLMARK pathway analysis for the BTNL9 high and low expression groups based on the TCGA database.
Estimating drug response for LUAD
Computational Analysis of Resistance (CARE) is a software that uses compound screening data to identify genome-scale biomarkers for targeted therapeutic strategies. Pearson correlation analysis between the gene expression profile of the cancer sample and the CARE scoring vector was used to group the patient as a responder or a non-responder .
Construction and evaluation of nomogram
We acquired TCGA LUAD RNA-seq data from the University of California, Santa Cruz (UCSC) Xena Browser (https://xenabrowser.net/). After screening, samples with missing clinical data and 0 days overall survival time were excluded, and a total of 501 samples were included. Next, we randomly divided the TCGA-LUAD cohort (n = 501) in a 7 to 3 ratio into a training (n = 352) and testing dataset (n = 149). We performed the R package “rms” to construct a nomogram based on the TNM stage and expression profile of BTNL9 using the training dataset. To evaluate the usefulness of the nomogram, the R package “ROCsurvival” was used to construct ROC for the prediction of the 1-, 3- and 5- year OS. The R package “ggDCA” was executed to create a decision analysis curve to evaluate the clinical utility of the nomogram. Finally, R package “rms” was applied to perform a calibration curve to evaluate the precision for predicting 1-, 3- and 5-year OS prediction of the LUAD cohort.
The relationship between BTNL9 expression and single cancer cell biological behavior of LUAD was determined using Pearson correlation analysis and Spearman’s correlation analysis of the correlation between BTNL9 and tumor mutation burden (TMB). In all the studies, P < 0.05 was considered statistically significant.
The high expression level of BTNL9 was associated with favorable survival of LUAD
Gene expression profile of BTNs, including BTN1A1, BTN2A1, BTN2A2, BTN2A3P, BTN3A1, BTN3A2, BTN3A3, BTNL2, BTNL3, BTNL8, BTNL9, BTNL10, and SKINTL was evaluated in normal and tumor lung tissues (Fig. 2A). Analysis showed that expression levels of BTNL8 and BTNL9 were significantly lower in tumor tissues compared with that of normal tissues (Fig. 2A). Furthermore, the expression level of BTNL9 was significantly negatively correlated with the clinical stage, lymph node metastasis stage, and p53 mutation. Concurrently, the expression level of BTNL8 was significantly negatively correlated with the clinical stage and N stage (Fig. 2B). However, survival analysis showed that BTNL8 was not significantly correlated with OS, whereas BTNL9 was significantly correlated with OS in LUAD (Fig. 2C). Validation of the prognosis value of BTNL9 in LUAD cohorts using PrognoScan  showed that BTNL9 expression was significantly correlated with RFS and OS in the GSE31210 dataset (n = 204). In addition, the expression level of BTNL9 was significantly associated with OS in GSE3141 cohort (n = 111) (Supplementary Table 1). Analysis using OncoLnc, UALCAN, and KM plotter showed that high expression of BTNL9 is significantly associated with better OS in LUAD (Fig. 2D). Although the two survival curves crossover occurred after 150 months (Fig. 2C), it is well beyond 5-years (60 months), and survival curves in the verification databases didn’t show crossover. Thus, we considered the results in this study reliable and stable. These findings imply that BTNL9 is a critical immune checkpoint of BTNs in LUAD.
Furthermore, the clinical distribution of BTNL9 was explored. TCGA-LUAD dataset was divided into the high and low groups based on the median gene expression level of BTNL9. The clinical characteristics, including age, gender, race, TNM staging, ECOG score, EGFR mutations, KRAS mutations, and radiotherapy information, were compared between the two groups. Analysis showed no significant difference in these clinical characteristics between the two groups (data not shown).
Pan- cancer gene expression and prognostic value of BTNL9
To further understand BTNL9 expression in pan-cancer, analysis of the dataset was performed using the Oncomine database. The findings showed that BTNL9 expression level was significantly lower in breast cancer, one colon cancer cohort, lung cancer, kidney cancer, and crabtree uterus cancer than normal tissues. However, BTNL9 expression was significantly higher in the brain and CNS cancer, colorectal cancer, esophageal cancer, leukemia, and lymphoma, than in normal tissue (Fig. 3A and Supplementary Table 2). Further, the BTNL9 expression profile was explored using TCGA RNA sequencing data (TIMER). The BTNL9 expression level was significantly downregulated in Bladder Urothelial Carcinoma (BLCA), Breast invasive carcinoma (BRCA), Cholangiocarcinoma (CHOL), Esophageal carcinoma (ESCA), Glioblastoma multiforme (GBM), Head and Neck squamous cell carcinoma (HNSC), Kidney renal papillary cell carcinoma (KIRP), LUAD, and LUSC compared with normal tissues. In contrast, BTNL9 was significantly increased in Colon adenocarcinoma (COAD), Kidney Chromophobe (KICH), and Kidney renal clear cell carcinoma (KIRC) compared with normal tissues (Fig. 3B). Analysis using the TissGDB database gave a correlation coefficient of BTNL9 expression with LUAD’s RFS HR of 0.87 [95% CI (0.79, 0.96)], and that for the correlation between BTNL9 expression and OS HR was 0.87 [95% CI (0.8,0.96)] (Fig. 3C, D).
Upstream and downstream regulatory network of BTNL9
A previous study reports that the expression level of BTNL9 in LUAD is significantly lower than that in normal tissues. DNA methylation is a biological process through which methyl groups are added to DNA molecules by Methyltransferases (DNMTs). DNA methylation of a gene promoter region functions by inhibiting gene transcription. Correlation analysis between BTNL9 expression and DNA methylation marker DNMTs (including DNMT1, DNMT2, DNMT3A, DNMT3B) was conducted using the GEPIA tool. Analysis showed that expression of BTNL9 in normal lung tissue was positively correlated with DNMTs (r = 0.35, P = 0.0059); however, there was no correlation with DNMTs in LUAD (r = − 0.019, P = 0.67) (Fig. 4A, B). These findings show that DNA methylation may be involved in the pathogenesis of LUAD. To further explore the upstream regulation mechanisms of the BTNL9 expression, miRNAs that bind to BTNL9 were predicted by using miRMap , TargetScan , and miRWalk  databases. A total of 248 miRNAs common predicted miRNAs from the three databases were obtained (Fig. 4C) and used starBase  to validate the predicted binding miRNAs. Analysis showed that, hsa-miR-30b-3p, hsa-miR-4709-3p and hsa-miR-6514-3p were significantly positively correlated with BTNL9 expression (r = 0.312, P = 5.25E-13, r = 0.277, P = 1.74E-10, and r = 0.103, P = 0.02, respectively, Fig. 4D). In addition, the three miRNAs were highly expressed and significantly correlated with higher OS of LUAD patients (HR = 0.66, P = 0.0058, HR = 0.63, P = 0.0023, and HR = 0.73, P = 0.036, respectively) (Fig. 4D). Although the two survival curves of hsa-miR-4709-3p and hsa-miR-6514-3p crossover occurred after 100 months, it is well beyond 5-years (60 months). Thus, we considered the results in this study reliable. Furthermore, 18 lncRNAs were predicted to bind to BTNL9 using LncMap  database (Supplementary Table 3). These findings were verified using LncACTdb2.0  database. LncRNA AP001462.6 was predicted to bind to BTNL9, and the high expression level of AP001462.6 was significantly correlated with a high OS of LUAD patients (P = 0.049) (Fig. 4E).
Moreover, proteins implicated in binding BTNL9 were analyzed using the STRING  database. Analysis showed a total of 7 proteins that bind to BTNL9 (Fig. 4F). The top 2 binding proteins, including HTRA4 and TM4SF19, were predicted using the cytoHubba module of Cytoscape  (Fig. 4G). HTRA4 gene encodes a member of the HtrA protease family. HTRA4 plays a role as a secreted oligomeric chaperone protease to degrade misfolded secretory proteins . We hypothesized that low expression of BTNL9 in LUAD might be related to degradation through ubiquitination. Analysis using Ubibrowser  database showed that E3 (MARCH8) ligases could bind the substrate BTNL9 (Supplementary Table 4). In addition, BTNL9 has a potential E3 recognizing domain and two potential E3 identifying motifs (Fig. 4H, I).
Low expression of BTNL9 significantly enriches proteasome and increases cancer malignancy
Gene Set Enrichment Analysis (GSEA) analysis for KEGG and HALLMARK was performed using the Sangerbox tool to explore the two groups’ biological pathways. The findings showed that the top 3 significantly enriched KEGG pathways in the high BTNL9 expression group were vascular smooth muscle contraction, phosphatidylinositol signaling system, and abc transporters (Fig. 5A). On the other hand, the top 4 significantly enriched KEGGs pathways in the low BTNL9 expression group were pathways implicated in Parkinson’s disease, oxidative phosphorylation, DNA replication, and proteasome pathways (Fig. 5B). GSEA for the HALLMARK pathway showed that the top 3 pathways associated with high BTNL9 expression were bile acid metabolism, heme metabolism, and Wnt/beta-catenin signaling pathways. Further, the top 4 pathways associated with low BTNL9 expression were E2F targets, glycolysis, myc targets v1, and mTORC1 signaling (Fig. 5C, D). These findings imply that BTNL9 is involved in LUAD metabolic reprogramming.
Metabolic reprogramming is a hallmark of cancer, and intrinsic and extrinsic factors contribute to various metabolic phenotypes in tumors. As cancer develops from pre-tumor lesions to local, clinically obvious malignant tumors to metastatic cancer, metabolism changes the phenotype and dependence . Single-cell RNA (scRNA) analysis of LUAD using CancerSEA  database (Fig. 5E) showed that BTNL9 expression is significantly negatively correlated with tumor malignant features including invasion (r = − 0.53, P < 0.0001), metastasis (r = − 0.35, P = 0.011), EMT (r = − 0.47, P = 0.0006), proliferation (r = − 0.37, P = 0.0086), Hypoxia (r = − 0.36, P = 0.011), and DNA damage (r = − 0.34, P = 0.017) (Fig. 5F). This finding implies that low expression of BTNL9 is significantly associated with the malignant features of LUAD.
Correlation between BTNL9 and infiltrating immune cell markers
Spearman’s correlation analysis of the correlation between expression of BTNL9 and tumor mutation burden (TMB) in the TCGA LUAD cohort showed that BTNL9 is significantly negatively correlated with TMB (P = 1.4E-9) (Fig. 6A). Analysis of somatic mutation pattern of BTNL9 in LUAD using the SangerBox tool showed that the mutation frequency of BTNL9 in LUAD was 1.41% (Fig. 6B). Genetic mutations are implicated in the tumor microenvironment (TME) ; therefore, the relationship between the expression of BTNL9 and the immune score was determined using the ESTIMATE algorithm in the SangerBox tool. Analysis showed that BTNL9 was significantly positively correlated with ImmuneScore (r = 0.129, P = 0.003) and ESTIMATEScore (r = 0.106, P = 0.016). However, the expression of BTNL9 was not significantly correlated with StromalScore (Fig. 6C-E).
Analysis of the correlation between gene expression of BTNL9 and infiltrating level of immune cells using TIMER  database showed that BTNL9 was negatively correlated with tumor purity (r = − 0.124, P = 5.6E-03). On the other hand, gene expression of BTNL9 was significantly positively correlated with B cells (r = 0.24, P = 8.88E-8), CD4+T (r = 0.283, P = 2.24E-10) and macrophages (r = 0.209, P = 3.35E-6) (Fig. 6F). Moreover, survival analysis showed that high expression of BTNL9 in B cells (P = 0.000) and DC cells (P = 0.048) was correlated with better OS for LUAD (Fig. 6G).
A detailed analysis of TME infiltrated DC and B cells using the TIMER database showed that DC and its subtypes cDCs1 and cDCs2  are associated with BTNL9 expression before and after purity adjustment. GEPIA database analysis showed that normal lung tissue was not correlated with DC and its subtypes cDCs1 and cDCs2. However, DC and its subtypes were significantly positively correlated with LUAD (Table 1), implying that DCs regulated by BTNL9 may participate in LUAD immune response. B cells are heterogeneous and include two subtypes: naïve B cells and plasma B cells . TIMER analysis showed that total B cells and naïve B cells were significantly correlated with BTNL9 expression before and after purity adjustment. However, plasma B cells were not associated with BTNL9 expression before and after purity adjustment. GEPIA analysis showed that total B cells and naïve B cells were not correlated with BTNL9 expression in normal lung tissues; however, they were significantly positively correlated with BTNL9 expression in LUAD tissues. Plasma B cells showed no correlation with BTNL9 in both normal tissues and LUAD tissues (Table 1), indicating that BTNL9 may play a role in promoting naïve B cell antitumor immune response.
High expression of BTNL9 is associated with tyrosine kinase inhibitors response
BTNL9 expression was significantly positively correlated with CARE scores for several compounds retrieved from Cancer Cell Line Encyclopedia (CCLE), Genomics of Drug Sensitivity in Cancer (GDSC, previously named CGP), and The Cancer Therapeutics Response Portal (CTRP) cohorts, mainly including antiangiogenic tyrosine kinase inhibitors Axitinib, Nilotinib, Sorafenib, Pazopanib, Masitinib, and Ki8751 (Fig. 7, and Table 2). These findings show that immune checkpoint inhibitors based on BTNL9 plus antiangiogenic tyrosine kinase inhibitors could be developed as a potential chemotherapy-free combination treatment strategy for LUAD.
Independent predictive power of BTNL9 based on multivariate analysis
We used the R package “survival V3.2–10” to construct a Cox model, including known important clinical variables for OS, such as TNM stage, primary therapy outcome, and BTNL9 expression. We also used multivariate analysis to explore whether BTNL9 expression was an independent OS factor for TCGA-LUAD patients. The results demonstrated that higher BTNL9 expression significantly (p = 0.049) and independently increased OS (HR = 0.67, 95% CI 0.45–0.99) (Table 3).
Development of a nomogram predicting OS
A nomogram predicting the 1-, 3- and 5- year OS for TCGA-LUAD was constructed based on BTNL9 expression and TNM stage (Fig. 8A). We built the ROC for the training dataset and the testing dataset and calculated the area under the ROC (AUC) to validate the accuracy of the nomogram. The AUCs for 1-, 3- and 5-year OS were 0.642, 0.645, and 0.607 in the training set (Fig. 8B); 0.727, 0.545, and 0.631 in the testing set (Fig. 8C). These results suggested that the nomogram showed a consistent accuracy in the training and testing dataset. We then conducted decision curve analysis (DCA) to evaluate the clinical usefulness, and the result showed that the nomogram provided an additional benefit compared to the “treat-all” and “treat-none” strategies in both the training and testing dataset (Fig. 8D-E). Finally, to compare the consistency of the model predictions with actual clinical outcomes, calibration curves for 1-year, 3-year, and 5-year OS were created for the training and testing dataset (supplementary Fig. 1A-F). The calibration curves showed consistent agreement between the predicted and observed values for 1-, 3- and 5-year OS.
Immune checkpoint inhibition or adoptive cell therapy has significantly changed the cancer treatment paradigm and has resulted in an age of chemotherapy-free NSCLC . Immune checkpoint inhibitors figure prominently in achieving chemotherapy-free cancer treatment. BTNs are immune checkpoints in several cancer types; however, the functions of BTNs have not been explored in LUAD. This study shows that BTNL9 is poorly expressed in LUAD tissues, and its low expression is correlated with a lower probability of 1, 3, 5-year OS based on a nomogram model. In addition, this team explored the mechanisms behind BTNL9 low expression. This study shows that mutated p53 results in a significant decrease in BTNL9 expression (Fig. 2B). Approximately 46% of LUAD patients possess p53 mutation . Breast cancer exhibits a low expression level of BTNL9, which can be targeted to inhibit proliferation and metastasis through the p53/CDC25C and p53/GADD45 signaling pathway . In addition, BTNL9 plays a role as a transcriptional modulator through epigenetic regulation and post-transcriptional modification.
DNA methylation is the most common form of DNA modification. It plays a vital role in normal cell physiology, and increased DNA methylation, and loss of demethylation, are observed in different cancer types. DNMTs are implicated in abnormal DNA methylation. Gene body hypermethylation activates oncogenes, and promotion of hypermethylation causes suppression of tumors . GEPIA analysis showed that BTNL9 and DNMTs correlate significantly with normal lung tissues but not LUAD tumorigenesis (Fig. 4A, B). However, the molecular mechanism of regulation of BTNL9 by DNMTs in LUAD has not been explored.
miRNA and lncRNA are non-coding RNAs involved in tumor promotion and suppression, depending on the tumor type . A total of 3 miRNAs (hsa-miR-30b-3p, hsa-miR-4709-3p, and hsa-miR-6514-3p) were significantly positively correlated BTNL9 in LUAD, and their high expression was significantly associated with longer OS. Previous studies reported that hsa-miR-30b-3p plays a role as an antitumor miRNA [42, 43]. Hsu Y-L et al. reported that BTNL9 acts as a tumor suppressor in LUAD and is regulated by hsa-miR-183-5p; however, the specific regulatory network was not reported . In addition, lncRNA AP001462.6 was shown to bind to BTNL9, and the high expression level of this lncRNA was significantly correlated with longer OS in LUAD patients (Fig. 4E).
Moreover, analysis of the protein interaction network of BTNL9 showed that the interacting proteins played a role in immune regulation, protease hydrolysis, and serine/threonine kinase regulation. Notably, protease hydrolysis is related to ubiquitination and degradation of proteins. Analysis showed that BTNL9 has a potential E3 recognizing domain binding site (Fig. 4H, I). These findings indicate that BTNL9 in LUAD may be regulated by DNA methylation and non-coding RNA. In addition, BTNL9 protein may be held by ubiquitination and degradation after translation.
We performed a GSEA analysis to explore the biological function of BTNL9 in LUAD. Functional analysis showed low expression levels of BTNL9 in energy metabolism (oxidative phosphorylation, glycolysis, myc targets v1 , and mTORC1 signaling ), DNA replication, and protease hydrolysis. Metabolic reprogramming triggers selective gene amplification and a large gene family, which drives cellular functions to promote cancer cell growth and proliferation . The above functions were subsequently verified from a single cell perspective. Findings showed that BTNL9 was significantly negatively correlated with cancer cell malignant behaviors such as proliferation, invasion, EMT, metastasis, and hypoxia. This result indicates that BTNL9 may play a role in LUAD tumor suppression.
TME is a potential predictor of response to an immune checkpoint inhibitor. Analysis of the relationship between BTNL9 and TME showed that the mutation frequency of BTNL9 in LUAD was about 1.14%, and BTNL9 was significantly negatively correlated with TMB. Furthermore, BTNL9 was significantly positively correlated with ImmuneScore and ESTIMATEScore. Previous studies report that high ImmuneScore and ESTIMATEScore are positively associated with a good prognosis of LUAD . This finding shows that BTNL9 plays an important role in TME immune regulation. Moreover, the correlation between BTNL9 and TILs showed that BTNL9 was significantly negatively correlated with tumor purity, and previous studies report that low tumor purity is associated with poor prognosis . Although BTNL9 was significantly correlated with B, CD4 + T, and macrophages, survival analysis showed that BTNL9 was only significantly correlated with B cells and DC cells.
DCs act pivotally in shaping innate and adaptive immune responses because they have a unique ability to initiate T-cell responses and promote their differentiation into effector lineages . B cells play antigen presentation, cytotoxicity, and antibody production functions, which are essential in adaptive immunity . TIMER and GEPIA database analysis showed that the expression level of BTNL9 was not correlated with levels of DCs (cDC1s and cDC2s) in normal adjacent tissues; however, BTNL9 was significantly associated with levels of DCs (cDC1s and cDC2s) in LUAD tissues (Table 1). cDC1s can migrate to tumor-draining lymph nodes, activate and attract T cells, secrete cytokines, and present antigens in TME, promoting local cytotoxic T cells . cDC2s present antigens to MHC II, activate CD4 + T cells, and effectively polarize TILs into anti-tumor T helper cell 1 (Th1) or Th17 phenotype . The BTNL9 expression level was not correlated with B cells (naïve B cells) in normal adjacent tissues; however, it was significantly associated with levels of B cells (naïve B cells) in LUAD tissues, except for plasma B cells. This finding implies that BTNL9 regulates the function of naïve B cells in TME. Previous studies report that naïve B cells are down-regulated in advanced NSCLC and are correlated with poor prognosis . Furthermore, CARE database analysis showed that BTNL9 expression is associated with effective antiangiogenic tyrosine kinase inhibitors response (Fig. 7, and Table 2). More data are being awaited to confirm this preliminary observation.
Notably, this study had a few limitations. Firstly, our findings are entirely based on public databases using bioinformatics analysis, and therefore further molecular biology experiments should be performed to verify these results. In addition, all findings presented here were developed using Database algorithms. Secondly, BTNL9 may not be detected by histopathology due to a lack of LUAD tissue samples. Further, OS analyses can’t be performed based on histopathology results. Finally, the study did not verify the role of BTNL9 in predicting immune responses in LUAD patients due to the lack of clinical cohorts of immunotherapy-treated LUAD patients. Taken together, results and conclusions based entirely on bioinformatics are informative and can lay down the foundations for more robust studies, but they do not replace experimentation.
In summary, the findings of this study show an association between immune checkpoint BTNL9 and OS in LUAD patients. Transcriptional regulation and post-transcriptional regulation are potential mechanisms for down-regulating BTNL9 expression, resulting in more malignant biological characteristics in LUAD. BTNL9 may modify the TME by enrichment of naïve B cells and DCs and promoting immune response and antiangiogenic tyrosine kinase inhibitors response.
Availability of data and materials
All data in this study were retrieved from public and open-source databases. All databases and its web address were listed below: GEPIA (http://gepia2.cancer-pku.cn/#index), Oncomine (https://www.oncomine.org/), KM plotter (http://www.kmplot.com), UALCAN (http://ualcan.path.uab.edu/index.html), OncoLnc (http://www.oncolnc.org/), TissGDB (https://bioinfo.uth.edu/TissGDB/index.html?csrt=17567836353851237153), TIMER (http://timer.comp-genomics.org/), Sangerbox (http://sangerbox.com/Index), miRMap (https://mirmap.ezlab.org/), TargetScan (http://www.targetscan.org/vert_72/), miRWalk (http://mirwalk.umm.uni-heidelberg.de/), starBase (http://starbase.sysu.edu.cn/panCancer.php), CancerSEA (biocc.hrbmu.edu.cn/CancerSEA/goSearch).
Lung squamous cell carcinoma
Non-small cell lung cancer
- BTN :
Gene expression profiling interactive analysis
Tumor immune estimation resource
Computational Analysis of REsistance
Gene expression omnibus
Eastern Cooperative Oncology Group
Gene set enrichment analysis
Single nucleotide mutations
Tumor mutation burden
Cancer Cell Line Encyclopedia
Genomics of Drug Sensitivity in Cancer
Cancer Therapeutics Response Portal
Tumor-infiltrating immune cells
Chen W, Zheng R, Baade PD, Zhang S, Zeng H, Bray F, et al. Cancer statistics in China, 2015. CA Cancer J Clin. 2016;66(2):115–32. https://doi.org/10.3322/caac.21338.
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–424. https://doi.org/10.3322/caac.21492.
Lu T, Yang X, Huang Y, Zhao M, Li M, Ma K, et al. Trends in the incidence, treatment, and survival of patients with lung cancer in the last four decades. Cancer Manag Res. 2019;11:943–53. https://doi.org/10.2147/CMAR.S187317.
Garon EB, Hellmann MD, Rizvi NA, Carcereny E, Leighl NB, Ahn MJ, et al. Five-year overall survival for patients with advanced non–small-cell lung Cancer treated with Pembrolizumab: results from the phase I KEYNOTE-001 study. J Clin Oncol. 2019;37(28):2518–27. https://doi.org/10.1200/JCO.19.00934.
Davern M, Lysaght J. Cooperation between chemotherapy and immunotherapy in gastroesophageal cancers. Cancer Lett. 2020;495:89–99. https://doi.org/10.1016/j.canlet.2020.09.014.
Guo Y, Wang AY. Novel Immune Check-Point Regulators in Tolerance Maintenance. Front Immunol. 2015;6:421. https://doi.org/10.3389/fimmu.2015.00421.
Abeler-Dörner L, Swamy M, Williams G, Hayday AC, Bas A. Butyrophilins: an emerging family of immune regulators. Trends Immunol. 2012;33(1):34–41. https://doi.org/10.1016/j.it.2011.09.007.
Malinowska M, Tokarz-Deptuła B, Deptuła W. Butyrophilins: an important new element of resistance. Cent-Eur J Immunol. 2017;42(4):399–403. https://doi.org/10.5114/ceji.2017.72806.
Arnett HA, Viney JL. Immune modulation by butyrophilins. Nat Rev Immunol. 2014;14(8):559–69. https://doi.org/10.1038/nri3715.
Zhou C. Lung cancer molecular epidemiology in China: recent trends. Transl Lung Cancer Res. 2014;3(5):270–9. https://doi.org/10.3978/j.issn.2218-6751.2014.09.01.
Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 2017;45(W1):W98–W102. https://doi.org/10.1093/nar/gkx247.
Clough E, Barrett T. The Gene Expression Omnibus Database. Methods Mol Biol. 2016;1418:93–110.
Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–w514. https://doi.org/10.1093/nar/gkaa407.
Győrffy B, Surowiak P, Budczies J, Lánczky A. Online survival analysis software to assess the prognostic value of biomarkers using transcriptomic data in non-small-cell lung cancer. PLoS One. 2013;8(12):e82241. https://doi.org/10.1371/journal.pone.0082241.
Chandrashekar DS, Bashel B, Balasubramanya SAH, Creighton CJ, Ponce-Rodriguez I, Chakravarthi B, et al. UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia. 2017;19(8):649–58.
Anaya J. OncoLnc: linking TCGA survival data to mRNAs, miRNAs, and lncRNAs. PeerJ Computer Science. 2016;2:e67. https://doi.org/10.7717/peerj-cs.67.
Rhodes DR, Yu J, Shanker K, Deshpande N, Varambally R, Ghosh D, et al. ONCOMINE: a cancer microarray database and integrated data-mining platform. Neoplasia. 2004;6(1):1–6.
Kim P, Park A, Han G, Sun H, Jia P, Zhao Z. TissGDB: tissue-specific gene database in cancer. Nucleic Acids Res. 2017;46(D1):D1031–8. https://doi.org/10.1093/nar/gkx850.
Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 2019;47(W1):W556–60. https://doi.org/10.1093/nar/gkz430.
Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 2017;77(21):e108–10. https://doi.org/10.1158/0008-5472.CAN-17-0307.
Dai D, Chen B, Feng Y, Wang W, Jiang Y, Huang H, et al. Prognostic value of prostaglandin I2 synthase and its correlation with tumor-infiltrating immune cells in lung cancer, ovarian cancer, and gastric cancer. Aging. 2020;12(10):9658–85. https://doi.org/10.18632/aging.103235.
Vejnar CE, Blum M, Zdobnov EM. miRmap web: Comprehensive microRNA target prediction online. Nucleic Acid Res. 2013;41(Web Server issue):W165–8.
Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. eLife. 2015;4:e05005. https://doi.org/10.7554/eLife.05005.
Sticht C, De La Torre C, Parveen A, Gretz N. miRWalk: an online resource for prediction of microRNA binding sites. PLoS One. 2018;13(10):e0206239-e0206239. https://doi.org/10.1371/journal.pone.0206239.
Li J-H, Liu S, Zhou H, Qu L-H, Yang J-H. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein–RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 2013;42(D1):D92–7. https://doi.org/10.1093/nar/gkt1248.
Jiang P, Lee W, Li X, Johnson C, Liu JS, Brown M, et al. Genome-Scale Signatures of Gene Interaction from Compound Screens Predict Clinical Efficacy of Targeted Cancer Therapies. Cell Syst. 2018;6(3):343–354.e345.
Mizuno H, Kitada K, Nakai K, Sarai A. PrognoScan: a new database for meta-analysis of the prognostic value of genes. BMC Med Genet. 2009;2(1):18. https://doi.org/10.1186/1755-8794-2-18.
Li Y, Li L, Wang Z, Pan T, Sahni N, Jin X, et al. LncMAP: Pan-cancer atlas of long noncoding RNA-mediated transcriptional network perturbations. Nucleic Acids Res. 2018;46(3):1113–23. https://doi.org/10.1093/nar/gkx1311.
Wang P, Li X, Gao Y, Guo Q, Wang Y, Fang Y, et al. LncACTdb 2.0: an updated database of experimentally supported ceRNA interactions curated from low-and high-throughput experiments. Nucleic Acids Res. 2019;47(D1):D121–7. https://doi.org/10.1093/nar/gky1144.
Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2018;47(D1):D607–13. https://doi.org/10.1093/nar/gky1131.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. https://doi.org/10.1101/gr.1239303.
Li Y, Xie P, Lu L, Wang J, Diao L, Liu Z, et al. An integrated bioinformatics platform for investigating the human E3 ubiquitin ligase-substrate interaction network. Nat Commun. 2017;8(1):347. https://doi.org/10.1038/s41467-017-00299-9.
Faubert B, Solmonson A, DeBerardinis RJ. Metabolic reprogramming and cancer progression. Science. 2020;368(6487):eaaw5473.
Yuan H, Yan M, Zhang G, Liu W, Deng C, Liao G, et al. CancerSEA: a cancer single-cell state atlas. Nucleic Acids Res. 2018;47(D1):D900–8. https://doi.org/10.1093/nar/gky939.
Giraldo NA, Sanchez-Salas R, Peske JD, Vano Y, Becht E, Petitprez F, et al. The clinical role of the TME in solid cancer. Br J Cancer. 2019;120(1):45–53. https://doi.org/10.1038/s41416-018-0327-z.
Brown CC, Gudjonson H, Pritykin Y, Deep D, Lavallée V-P, Mendoza A, et al. Transcriptional Basis of Mouse and Human Dendritic Cell Heterogeneity. Cell. 2019;179(4):846–863.e824.
Chen J, Tan Y, Sun F, Hou L, Zhang C, Ge T, et al. Single-cell transcriptome and antigen-immunoglobin analysis reveals the diversity of B cells in non-small cell lung cancer. Genome Biol. 2020;21(1):152. https://doi.org/10.1186/s13059-020-02064-6.
Chu T, Zhong R, Zhong H, Zhang B, Zhang W, Shi C, et al. Phase Ib study of Sintilimab plus Anlotinib as first-line therapy in patients with advanced non-small cell lung Cancer. J Thorac Oncol. 2021;16(4):643-52. https://doi.org/10.1016/j.jtho.2020.11.026.
Collisson EA, Campbell JD, Brooks AN, Berger AH, Lee W, Chmielecki J, et al. Comprehensive molecular profiling of lung adenocarcinoma. Nature. 2014;511(7511):543–50. https://doi.org/10.1038/nature13385.
Mo Q, Xu K, Luo C, Zhang Q, Wang L, Ren G. BTNL9 is frequently downregulated and inhibits proliferation and metastasis via the P53/CDC25C and P53/GADD45 pathways in breast cancer. Biochem Biophys Res Commun. 2021;553:17–24. https://doi.org/10.1016/j.bbrc.2021.03.022.
Cheng Y, He C, Wang M, Ma X, Mo F, Yang S, et al. Targeting epigenetic regulators for cancer therapy: mechanisms and advances in clinical trials. Signal Transduct Target Ther. 2019;4(1):62. https://doi.org/10.1038/s41392-019-0095-0.
Zhong K, Chen K, Han L, Li B. MicroRNA-30b/c inhibits non-small cell lung cancer cell proliferation by targeting Rab18. BMC Cancer. 2014;14(1):703. https://doi.org/10.1186/1471-2407-14-703.
Gao D, Zhou Z, Huang H. miR-30b-3p Inhibits Proliferation and Invasion of Hepatocellular Carcinoma Cells via Suppressing PI3K/Akt Pathway. Front Genet. 2019;10:1274. https://doi.org/10.3389/fgene.2019.01274.
Hsu Y-L, Hung J-Y, Lee Y-L, Chen F-W, Chang K-F, Chang W-A, et al. Identification of novel gene expression signature in lung adenocarcinoma by using next-generation sequencing data and bioinformatics analysis. Oncotarget. 2017;8(62):104831–54. https://doi.org/10.18632/oncotarget.21022.
Miller DM, Thomas SD, Islam A, Muench D, Sedoris K. C-Myc and cancer metabolism. Clin Cancer Res. 2012;18(20):5546–53. https://doi.org/10.1158/1078-0432.CCR-12-0977.
Valvezan AJ, Manning BD. Molecular logic of mTORC1 signalling as a metabolic rheostat. Nat Metab. 2019;1(3):321–33. https://doi.org/10.1038/s42255-019-0038-7.
Öjlert ÅK, Halvorsen AR, Nebdal D, Lund-Iversen M, Solberg S, Brustugun OT, et al. The immune microenvironment in non-small cell lung cancer is predictive of prognosis after surgery. Mol Oncol. 2019;13(5):1166–79. https://doi.org/10.1002/1878-0261.12475.
Mao Y, Feng Q, Zheng P, Yang L, Liu T, Xu Y, et al. Low tumor purity is associated with poor prognosis, heavy mutation burden, and intense immune phenotype in colon cancer. Cancer Manag Res. 2018;10:3569–77. https://doi.org/10.2147/CMAR.S171855.
Hu X, Zhang J, Wang J, Fu J, Li T, Zheng X, et al. Landscape of B cell immunity and related immune evasion in human cancers. Nat Genet. 2019;51(3):560–7. https://doi.org/10.1038/s41588-018-0339-x.
Böttcher JP, Reis e Sousa C. The role of type 1 conventional dendritic cells in Cancer immunity. Trends Cancer. 2018;4(11):784–92. https://doi.org/10.1016/j.trecan.2018.09.001.
Perez CR, De Palma M. Engineering dendritic cell vaccines to improve cancer immunotherapy. Nat Commun. 2019;10(1):5408–5408. https://doi.org/10.1038/s41467-019-13368-y.
The authors thank the developers who provided access for the GEPIA, TIMER, KM plotter, UALCAN, OncoLnc, Oncomine, TissGDB, starBase, CancerSEA, and Sangerbox database.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The original version of this article was revised: Figure 2 has been update. The gene BTN3A1 had Chinese text characters, and BTN3A2, BTN3A3, BTNL2, BTNL3, BTNL8, BTNL9, BTNL10, and SKINTL were not marked.
The calibration curves for predicting 1-year, 3-year, and 5-year OS in LUAD. (A, B, and C) Calibration curves for 1, 3, 5-year OS in training dataset; (D, E, and F) Calibration curves for 1, 3, 5-year OS in testing dataset.
About this article
Cite this article
Ma, W., Liang, J., Mo, J. et al. Butyrophilin-like 9 expression is associated with outcome in lung adenocarcinoma. BMC Cancer 21, 1096 (2021). https://doi.org/10.1186/s12885-021-08790-9
- Butyrophilin-like 9
- B cells
- Dendritic cells
- Lung adenocarcinoma