Butyrophilin-like 9 expression is associated with outcome in lung adenocarcinoma

Background 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. Methods 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. Results 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. Conclusions 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. Supplementary Information The online version contains supplementary material available at 10.1186/s12885-021-08790-9.


Background
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%) [3]. 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% [4]. Immune checkpoint inhibitors (ICIs) block immune checkpoint signaling, thus alleviating antitumor immunity, and significantly improving five year-OS in NSCLC.
Lung adenocarcinoma (LUAD) has been the most prevalent histopathological subtype of NSCLC in China since 2014 [10]. In this study, we explored the relationship between the expression level of BTNs and LUAD prognosis. Significant survival-related BTNs were screened using GEPIA [11]. Datasets used for analysis in this study were retrieved from Gene Expression Omnibus [12], TIMER [13], KM plotter [14], UALCAN [15], OncoLnc [16], Oncomine [17], TissGDB [18] 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 [19]. 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 [17] using default settings with a P-value of 0.01, a fold change of 1.5, and a top 10% gene ranking.

Survival analysis
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 [14], UALCAN [15], and OncoLnc databases [16]. TissGDB is a tissue-specific gene annotation database in Fig. 1 A flow chart of the study design cancer [18]. 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 [21]. 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 [22], TargetScan [23], and miRWalk [24] were used to predict miRNAs that can bind to BTNL9. Predicted miRNAs obtained from the three databases were further verified using the starBase database [25].
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 [26].

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.

Statistical analysis
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.

Results
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 [27] 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  Correlation between BTNL9 expression and methyltransferases (DNMTs) such as DNMT1, DNMT2, DNMT3A, and DNMT3B in LUAD and adjacent tissues using GEPIA. (C) Predicted miRNAs that bind to BTNL9 using miRMap, TargetScan, and miRWalk databases presented as a Venn diagram. (D) Overlapping 248 miRNAs verified using the StarBase database, hsa-miR-30b-3p, hsa-miR-4709-3p, and hsa-miR-6514-3p were screened. (E) Predicted LncRNAs that bind to BTNL9 were predicted using the LncMap database. AP001462.6 was verified and screened using the LncACTdb2.0 database. (F, G) BTNL9 interacting proteins were identified using the STRING database and edited and visualized using Cytoscape software (V3.7.2). Hub genes were screened using the cytoHubba module in Cytoscape. (H, I) Ubibrowser database predicts that the substrate BTNL9 can be bound by E3 (MARCH8) ligases, with one potential E3 recognizing domain and two potential E3 identifying motifs 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, miR-NAs that bind to BTNL9 were predicted by using miR-Map [22], TargetScan [23], and miRWalk [24] databases. A total of 248 miRNAs common predicted miRNAs from the three databases were obtained (Fig. 4C) and used starBase [25] to validate the predicted binding miR-NAs. 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,  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 [28] database (Supplementary Table 3). These findings were verified using LncACTdb2.0 [29] 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 [30] 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 [31] (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 [19]. We hypothesized that low expression of BTNL9 in LUAD might be related to degradation through ubiquitination. Analysis using Ubibrowser [32] 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 [33]. Single-cell RNA (scRNA) analysis of LUAD using CancerSEA [34] (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) [35]; 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 [13] 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 [36] 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 [37]. 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, 3year, 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.

Discussion
Immune checkpoint inhibition or adoptive cell therapy has significantly changed the cancer treatment paradigm  and has resulted in an age of chemotherapy-free NSCLC [38]. 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 [39]. 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 [40]. 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 [41]. 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 [41]. 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 [44]. 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 [45], and mTORC1 signaling [46]), 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 [45]. 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 [47]. 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 [48]. 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 [36]. B cells play antigen presentation, cytotoxicity, and antibody production functions, which are essential in adaptive immunity [49]. 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 [50]. 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 [51]. 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 [37]. 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.

Conclusion
In summary, the findings of this study show an association between immune checkpoint BTNL9 and OS in LUAD patients. Transcriptional regulation and posttranscriptional 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.