BTNL9 mRNA expression in LUAD
The "Tidyverse" and "Table 1" packages of R software (version 4.0.1) was used to extract the clinical parameters of the LUAD data set in TCGA. Further, we analyzed the age, gender, race, TNM staging, ECOG score, EGFR mutations, KRAS mutations, and radiotherapy information of the two groups of patients with high and low expression of BTNL9. There is no significant difference in two groups (data not shown), and we organized them into a table for descriptive analysis (Table 1).
Btnl9 Is The Key Gene Of Butyrophilins In Luad
The expression of BTNL8 and BTNL9 were significantly different in LUAD, and both tumor tissues were significantly lower than normal tissues (Fig. 1A). BTNL9 significantly correlated with the clinical stage, N stage and P53 mutation, while BTNL8 significantly correlated with the clinical stage and N stage (Fig. 1B). However, survival analysis revealed that BTNL8 was not significantly correlated with OS, while BTNL9 was significantly correlated with OS in LUAD (Fig. 1C). PrognoScan[25] analysis of the NSCLC dataset found that BTNL9 expression was significantly correlated with RFS and OS in the GSE31210 study, and BTNL9 expression was significantly correlated with OS in GSE3141 study (Supplementary Table 1). KM plotter, UALCAN, and OncoLnc were used to validate the findings, and the results showed that high expression of BTNL9 significantly improved the OS of LUAD (Fig. 1D). Therefore, we regarded BTNL9 as the key gene of BTN and BTNL family in LUAD.
Expression Of Btnl9 In Pan-cancers And Its Prognostic Correlation
We analyzed the expression of BTNL9 in pan-cancer, and analysis of pan-cancers from the Oncomine database showed that BTNL9 expression was significantly reduced in breast cancer, one dataset of colon cancer, lung cancer, kidney cancer, and crabtree uterus cancer compared with the normal tissue. However, BTNL9 expression significantly increased in brain and CNS cancer, four datasets of colorectal cancer, esophageal cancer, leukemia, and lymphoma (Fig. 2A and Supplementary Table 2). We further examined BTNL9 expression using TCGA RNA sequencing data (TIMER). As shown in Fig. 2B, BTNL9 expression was significantly reducedin Bladder Urothelial Carcinoma (BLCA), Breast invasive carcinoma (BRCA), Cholangio carcinoma (CHOL), Esophageal carcinoma (ESCA), Glioblastoma multiforme (GBM), Head and Neck squamous cell carcinoma (HNSC), Kidney renal papillary cell carcinoma (KIRP2), LUAD, and LUSC, but significantly increasedin Colon adenocarcinoma (COAD), Kidney Chromophobe (KICH), and Kidney renal clear cell carcinoma (KIRC). The TissGDB database analysis revealed that the correlation of BTNL9 expression with LUAD's RFS HR was 0.87[ 95%CI (0.79, 0.96)], and the OS HR was 0.87 [95%CI (0.8,0.96)] (Fig. 2C, 2D).
Upstream And Downstream Regulatory Network Of Btnl9
In previous study, we found that the expression of BTNL9 in LUAD was significantly lower compared with the normal tissues. DNA methylation is related to gene expression. Methyltransferases (DNMTs) including DNMT1, DNMT2, DNMT3A, and DNMT3B. Analysis of the GEPIA database revealed that the expression of BTNL9 in normal lung tissue was significantly positively correlated with DNMTs (r = 0.35,P = 0.059), but there was no correlation in LUAD (r=-0.019,P = 0.67) (Fig. 3A, 3B), suggesting that DNA methylation may be involved in the pathogenesis of LUAD. We further analyzed the upstream regulation of BTNL9, predicted miRNAs that bind to BTNL9 through miRMap[21], TargetScan[26], and miRWalk[23], and considered the intersection of the 3 databases to obtain 248 miRNAs (Fig. 3C). The miRNAs were verified in Starbase[24] database and found that, hsa-miR-30b-3p, hsa-miR- 4709-3p and hsa-miR-6514-39 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), and the 3 miRNAs were highly expressed and significantly prolonged the OS of LUAD (HR = 0.66, P = 0.0058, HR = 0.63, P = 0.0023, and HR = 0.73, P = 0.036, respectively) (Fig. 3D). Furthermore, we predicted the lncRNA binding to BTNL9 through LncMap[27] database and obtained 18 LncRNAs (Supplementary Table 3), which were verified and predicted the OS correlation in these LncRNAs through LncACTdb2.0[28] database. LncRNA AP001462.6 was found to bind to BTNL9, and the high expression of AP001462.6 significantly prolonged the OS in LUAD (P = 0.049) (Fig. 3E).
Moreover, we analyzed the protein network that binds to BTNL9 through the STRING[29] database and found that there were 7 binding proteins (Fig. 3F). The top 2 binding proteins were predicted using the cytoHubba module of Cytoscape[30] (Fig. 3G). The proteins that bind to BTNL9 were predicted to be HTRA4 and TM4SF19. HTRA4 gene encodes a member of the HtrA family of proteases, which may function as a secreted oligomeric chaperone protease to degrade misfolded secretory proteins[18]. We hypothesized that the low expression of BTNL9 in LUAD may be related to ubiquitination degradation. Analysis ofthe Ubibrowser[31] database revealed that E3 (MARCH8) ligases can bind the substrate BTNL9 (Confidence level is high, and score = 0.805), and there is one potential E3 recognizing domain and two potential E3 recognizing motif (Fig. 3H, 3I).
Low BTNL9 expression significantly enriches proteasome and increases cancer malignancy
To determine the biological function of BTNL9 expression on LUAD, the Sangerbox tool was used to divide the TCGA samples into high and low BTNL9 expression groups. Gene set enrichment analysis (GSEA) for KEGG and HALLMARK pathways were performed. The results revealed that the top 3 KEGGs pathways with high expression of BTNL9 were significantly enriched in vascular_smooth_muscle_contraction, phosphatidylinositol_signaling_system, and abc_transporters (Fig. 4A). The top 4 KEGGs pathways associated with low expression of BTNL9 were parkinsons_disease, oxidative_phosphorylation, DNA replication, and proteasome (Fig. 4B). GSEA for the HALLMARK pathway revealed that the top 3 pathways associated with high BTNL9 expression were bile_acid_metabolism, heme_metabolism, and wnt_beta_catenin_signaling, while the top 4 pathways associated with low BTNL9 expression were E2F_targets, glycolysis, myc_ targets v1, and mTORC1_signaling (Fig. 4C, 4D).
Metabolic reprogramming is a hallmark of cancer, and intrinsic and extrinsic factors contribute to metabolic phenotypes in tumors. As cancer develops from pre-tumor lesions to local, clinically obvious malignant tumors to metastatic cancer, the phenotype and dependence of metabolism also develop[32]. Through single-cell RNA (scRNA) analysis of LUAD in CancerSEA[33] database (Fig. 4E), BTNL9 expression was found to be 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. 4F). Therefore,the low expression of BTNL9 promotes the malignant features of LUAD.
Correlation between BTNL9 and markers of infiltrating immune cells and its prognosis
Tumor mutation burden (TMB) is measured by the number of somatic mutations that occur at an average of 1 Mb in the coding region (exon region) of the tumor cell genome (non-synonymous mutations). The total number of synonymous mutations indicates that the mutation types mainly include single nucleotide mutations (SNV) and small fragments of insertion/deletion (Indel) and other forms of mutations. Spearman’s correlation test analyzed the correlation between BTNL9 and TMB in TCGA-LUAD and found that BTNL9 was significantly negatively correlated with TMB (P = 1.4E-9) (Fig. 5A). Further, analysis of the somatic mutation pattern of BTNL9 in LUAD (Frame_Shift_Del, Missense, and Splice), revealed that the mutation frequency of BTNL9 in LUAD was 1.14% (Fig. 5B). Genetic mutations shape TME[34]. SangerBox tool used the R software package ESTIMATE to analyze the correlation between BTNL9 and stromal score. The results showed that BTNL9 was significantly positively correlated with ImmuneScore (r = 0.129, P = 0.003) and ESTIMATEScore (r = 0.106, P = 0.016). However, StromalScore was not significantly different (Fig. 5C, 5D,5E).
The correlation between BTNL9 and TME infiltrating immune cells was analyzed using TIMER[10] database and found that BTNL9 was negatively correlated with tumor purity (r=-0.124, P = 5.6E-03), and 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. 5F). Moreover, survival analysis showed that high expression of BTNL9 in B cells (P = 0) and DC cells (P = 0.048) significantly prolonged the OS of LUAD(Fig. 5G).
DC also have a trend of difference, and high DC expression of BTNL9 significantly prolonged OS in LUAD (Fig. 5F, 5G). We conducted a detailed analysis of TME infiltrated DC and B cells using the TIMER database, which revealed that DC and its subtypes cDCs1 and cDCs2[35] were all associated with BTNL9 expression before and after purity adjustment. GEPIA database analysis found that normal lung tissue was not correlated with DC and its subtypes cDCs1 and cDCs2, but was significantly positively correlated with LUAD (Table 2), indicating that all DCs regulated by BTNL9 may participate in the LUAD immune response. B cells are heterogeneous, and include two subtypes: naïve B cells and plasma B cells[36]. TIMER analysis found that total B cells and naïve B cells were significantly related to BTNL9 expression before and after purity adjustment, but plasma B cells were not correlated to BTNL9 expression before and after purity adjustment. GEPIA analysis found that total B cells and naïve B cells were not correlated to BTNL9 expression in normal lung tissues, but significantly positively correlated in LUAD. Plasma B cells showed no correlation with BTNL9 in both normal tissues and LUAD (Table 2), indicating that BTNL9 may promote naïve B cell anti-tumor immune response.
Table 2
Correlation analysis between BTNL9 and relate gene set markers of significant innate and adaptive immunity cells in TIMER and GEPIA database
| | | TIMER | | GEPIA |
| Immune cell | Marker | None | purity | | Normal | Cancer |
| | | Spearman's ρ | P Value | Spearman's ρ | P Value | | Spearman correlation coefficient | P Value | Spearman correlation coefficient | P Value |
| DC | CD1C | 0.36 | 1.65E-17 | 0.34 | 1.03E-14 | | -0.022 | 0.87 | 0.33 | 1.5E-13 |
| | HLA-DPA1 | 0.23 | 1.64E-07 | 0.20 | 8.47E-06 | | | | | |
| | HLA-DPB1 | 0.29 | 8.39E-12 | 0.27 | 8.03E-10 | | | | | |
| | HLA-DQB1 | 0.21 | 1.31E-06 | 0.18 | 8.16E-05 | | | | | |
| | HLA-DRA | 0.18 | 4.54E-05 | 0.14 | 0.0020 | | | | | |
| | ITGAX | 0.24 | 1.63E-08 | 0.22 | 1.10E-06 | | | | | |
| | NCR1 | 0.04 | 0.35980601 | 0.01 | 0.8450 | | | | | |
| | NRP1 | 0.03 | 0.41981037 | 0.03 | 0.4973 | | | | | |
| cDC1s | CD8A | 0.03 | 0.55316497 | -0.03 | 0.5613 | | -0.028 | 0.83 | 0.22 | 1.20E-06 |
| | CLEC9A | 0.40 | 1.00E-21 | 0.38 | 1.16E-18 | | | | | |
| | XCR1 | 0.41 | 2.70E-22 | 0.39 | 7.37E-20 | | | | | |
| cDC2s | CLEC12A | 0.22 | 3.50E-07 | 0.18 | 3.61E-05 | | 0.18 | 0.16 | 0.45 | 7.1E-25 |
| | ESAM | 0.51 | 1.01E-35 | 0.50 | 2.89E-33 | | | | | |
| B cell | CD19 | 0.23 | 1.93E-07 | 0.20 | 5.30E-06 | | 0.017 | 9.00E-01 | 0.37 | 8.30E-17 |
| | FCER2 | 0.40 | 7.97E-21 | 0.38 | 2.18E-18 | | | | | |
| | MS4A1 | 0.34 | 6.45E-16 | 0.33 | 2.59E-14 | | | | | |
| | SDC1 | 0.07 | 0.0917 | 0.09 | 0.0523 | | | | | |
| Naïve B cell | CD19 | 0.23 | 1.93E-07 | 0.20 | 5.30E-06 | | 0.082 | 0.54 | 0.43 | 8.60E-23 |
| | CD22 | 0.50 | 2.44E-34 | 0.51 | 1.10E-34 | | | | | |
| | CD83 | 0.33 | 6.56E-15 | 0.31 | 2.09E-12 | | | | | |
| | MS4A1 | 0.34 | 6.45E-16 | 0.33 | 2.59E-14 | | | | | |
| | TCL1A | 0.25 | 1.57E-08 | 0.21 | 1.31E-06 | | | | | |
| Plasma B cell | CD38 | 0.04 | 0.3699 | 0.08 | 0.0746 | | -0.19 | 0.15 | 0.059 | 0.19 |
| | TNFRSF17 | 0.04 | 0.4053 | 0.08 | 0.0814 | | | | | |
Btnl9 Expression Is To Be Associated With Drug Response
Computational Analysis of Resistance (CARE) is a software that uses compound screening data to identify genome-scale biomarkers of targeted therapeutic response. According to the Pearson correlation between the gene expression profile of the cancer sample and the CARE scoring vector, the patient is predicted as a responder or a non-responder[37]. BTNL9 expression has significant positive CARE scores for many compounds in Cancer Cell Line Encyclopedia (CCLE), Genomics of Drug Sensitivity in Cancer (GDSC, previously named CGP), and The Cancer Therapeutics Response Portal (CTRP) cohorts, suggesting that loss of BTNL9 expression may promote drug resistance toward many targeted therapies (Fig. 6, and Table 3).
Table 3
Loss of BTNL9 expression may promote drug resistance toward many targeted therapies in CGP, CCLE, and CTRP cohorts
| Drug | Target | t-value | p-value |
| IOX2 | EGLN1 | 5.91996 | 4.56E-09 |
| PHA-793887 | CDK9 | 4.92172 | 1.02E-06 |
| OSI-027 | MTOR | 4.55428 | 5.99E-06 |
| Ispinesib | KIF11 | 4.42265 | 1.09E-05 |
| Nilotinib | ABL1 | 4.3886 | 1.30E-05 |
| Axitinib | PDGFRA | 4.33456 | 1.64E-05 |
| NG25 | MAP4K2 | 4.29746 | 1.92E-05 |
| Nilotinib | KIT | 4.26379 | 2.26E-05 |
| BMS345541 | IKBKB | 4.12048 | 4.13E-05 |
| GSK525762A | BRD2 | 3.92075 | 9.50E-05 |
| CAY10603 | HDAC6 | 3.91591 | 9.69E-05 |
| TubastatinA | HDAC6 | 3.89552 | 0.000105 |
| GSK525762A | BRD4 | 3.83187 | 0.000136 |
CGP dataset | PHA-793887 | CDK1 | 3.4184 | 0.000658 |
| Fluorouracil | TYMS | 3.09815 | 0.002008 |
| TPCA-1 | IKBKB | 3.09042 | 0.002061 |
| 1256580-46-7 | ALK | 3.01459 | 0.002646 |
| CAL-101 | PIK3CD | 2.99145 | 0.002852 |
| Belinostat | HDAC6 | 2.82403 | 0.004852 |
| AT7519 | CDK9 | 2.81016 | 0.00506 |
| CP-466722 | ATM | 2.78047 | 0.005542 |
| Enzastaurin | PRKCB | 2.76951 | 0.00573 |
| SB590885 | BRAF_V600E.Mutation | 2.6446 | 0.008337 |
| Vorinostat | HDAC6 | 2.54112 | 0.011233 |
| Nutlin-3 | MDM2 | -2.94731 | 0.003295 |
| Dasatinib | EPHA2 | -3.23912 | 0.001305 |
| Quizartinib | FLT3 | -3.33865 | 0.000877 |
| 1173900-33-8 | PIK3CB | -3.55986 | 0.00039 |
| Linifanib | FLT3 | -5.1822 | 2.71E-07 |
| 870483-87-7 | CSF1R | -9.60837 | 7.14E-21 |
| Panobinostat | HDAC1 | 3.29253 | 0.001066 |
| abraxane | TUBB | 2.83794 | 0.00473 |
CCLE dataset | Palbociclib | RB1 | 2.64989 | 0.008358 |
| Topotecan | TOP1 | 2.45344 | 0.014498 |
| Sorafenib | FLT3 | -4.25026 | 2.56E-05 |
| 9-Fluoro-11,17,21-trihydroxy-16-methylpregna-1,4-diene-3,20-dione | NR3C1 | 5.30047 | 1.51E-07 |
| abraxane | TUBB | 4.50938 | 7.51E-06 |
| Alisertib | AURKB | 4.50069 | 7.83E-06 |
| PAC-1 | CASP3 | 4.27611 | 2.16E-05 |
| 660868-91-7 | PLK1 | 4.21764 | 2.79E-05 |
| Gossypol | BCL2 | 4.19093 | 3.12E-05 |
| Decitabine | DNMT1 | 4.18042 | 3.24E-05 |
| 722544-51-6 | AURKB | 4.14645 | 3.74E-05 |
| Etoposide | TOP2B | 4.0066 | 6.76E-05 |
| 3,5-bis(4-methylbenzylidene)piperidin-4-one | USP13 | 3.92444 | 9.47E-05 |
| 180002-83-9 | CNR2 | 3.91953 | 9.68E-05 |
| CICLOPIROX | RRM1 | 3.90668 | 0.000102 |
| Dabrafenib | BRAF_V600E.Mutation | 3.84863 | 0.00014 |
| BRD-K62801835-001-01-0 | EZH2 | 3.80929 | 0.000151 |
| BI2536 | PLK1 | 3.69725 | 0.000233 |
| Nilotinib | ABL1 | 3.67464 | 0.000255 |
| Cerulenin | HMGCS1 | 3.66163 | 0.000268 |
| TW-37 | BCL2 | 3.585 | 0.000358 |
| zebularine | DNMT1 | 3.57912 | 0.000366 |
| Pevonedistat | NAE1 | 3.55057 | 0.000409 |
| BAS02002358 | GPER1 | 3.52883 | 0.000442 |
| KU-60019 | ATM | 3.52842 | 0.000443 |
| narciclasine | RHOA | 3.47702 | 0.000536 |
| Ki8751 | PDGFRA | 3.43765 | 0.000619 |
| 4ly1 | HDAC1 | 3.39584 | 0.000719 |
CTRP dataset | SCHEMBL12182311 | EIF4E | 3.35831 | 0.000822 |
| Belinostat | HDAC1 | 3.35122 | 0.000886 |
| Masitinib | PDGFRA | 3.32101 | 0.000941 |
| BIBR1532 | TERT | 3.31812 | 0.000949 |
| UNII-UZ77T1VFBM | BIRC5 | 3.30878 | 0.000985 |
| Nutlin-3 | MDM2 | 3.2982 | 0.001018 |
| CHEMBL2058177 | EIF4E | 3.29052 | 0.001045 |
| NSC373989 | MDM2 | 3.24218 | 0.001239 |
| Imatinib | ABL1 | 3.22213 | 0.001327 |
| Tacrolimus | PPP3CB | 3.13945 | 0.001761 |
| Olaparib | PARP1 | 3.07789 | 0.002159 |
| SMR001317659 | PDE4A | 3.03729 | 0.002467 |
| Axitinib | PDGFRA | 3.02572 | 0.002563 |
| Pubchem_92131101 | KIF11 | 3.01909 | 0.002617 |
| abraxane | TUBB1 | 3.01218 | 0.002678 |
| GSK461364 | PLK1 | 2.97879 | 0.002987 |
| Sorafenib | PDGFRA | 2.93438 | 0.003443 |
| Nutlin-3 | TP53 | 2.88952 | 0.003968 |
| SMR000068650 | S1PR2 | 2.88871 | 0.003992 |
| MK-1775 | WEE1 | 2.88827 | 0.003985 |
| BRD-K53855319-001-01-2 | SIRT1 | 2.87063 | 0.004217 |
| Apicidin | HDAC1 | 2.7849 | 0.005485 |
| TelomeraseInhibitorIX | TERT | 2.7825 | 0.005527 |
| Pluripotin | MAPK1 | 2.76988 | 0.005741 |
| PHA-793887 | CDK1 | 2.75543 | 0.005999 |
| Vorinostat | HDAC1 | 2.74595 | 0.006174 |
| I-BET151 | BRD4 | 2.73756 | 0.006332 |
| SCHEMBL12182311 | EIF4A2 | 2.68306 | 0.007448 |
| 112522-64-2 | HDAC2 | 2.66828 | 0.007855 |
| RG108 | DNMT1 | 2.65008 | 0.008212 |
| PF184 | IKBKB | 2.63052 | 0.008701 |
| Pazopanib | PDGFRA | 2.63026 | 0.008704 |
| JQ-1 | BRD4 | 2.59578 | 0.009615 |
| SCHEMBL13833318 | HDAC1 | 2.56094 | 0.010633 |
| CAY10603 | HDAC6 | 2.55261 | 0.010887 |
| GSK525762A | BRD4 | 2.51599 | 0.012071 |
| 4CA-0620 | PLK1 | 2.43094 | 0.01529 |
| BCP9000801 | MDM2 | 2.41963 | 0.015783 |
| Belinostat | HDAC6 | 2.37335 | 0.018128 |
| Tipifarnib | FNTA | 2.36404 | 0.018536 |
| prima-1 | TP53 | 2.35318 | 0.018879 |
| BMS345541 | IKBKB | 2.35317 | 0.018872 |
| Gemcitabine | RRM1 | 2.3315 | 0.019984 |
| SMR001317659 | PDE4B | 2.28061 | 0.022841 |
| Doxorubicin | TOP2B | 2.26369 | 0.023866 |
| Rigosertib | PIK3CA | 2.22431 | 0.026433 |
| 3,5-di-tert-butylchalcone | RARA | -2.40858 | 0.016249 |
| Dasatinib | EPHA2 | -2.95126 | 0.003261 |