Butyrophilin-like 9 regulates immune inltration and serves as a prognostic marker in lung adenocarcinoma

Background: Butyrophilin (BTN) and butyrophilin-like (BTLN) belong to immunoglobulin superfamily, and also pertain the B7 co-stimulatory molecules family, which has multiple roles in immune modulation. Whether the BTNL9 expression in lung adenocarcinoma (LUAD) correlates with outcome has not been evaluated. Materials and methods: Oncomine and GEPIA were used to analyze the BTNL9, while its mRNA expression in LUAD and corresponding adjacent tissues was investigated using TIMER. The clinical prognosis of BTNL9 was assessed in the GEPIA, Kaplan-Meier plotter, and OncoLnc. Besides, the correlation between BTNL9 and tumor-inltrating immune cells (TILs) was analyzed using TIMER and GEPIA. The correlation between BTNL9 expression and drug response was analyzed using CARE. Results: BNL9 expression was signicantly low in LUAD. Low BTNL9 expression was associated with poor OS, and its expression was found to be regulated by both epigenetic regulation and post-transcriptional modication. BTNL9 expression was signicant positively correlated with ImmuneScore and ESTIMATEScore. Moreover, BTNL9 expression was positively correlated with inltrating levels of B cells, CD4+T, and macrophages, but Kaplan-Meier analysis showed that BTNL9 expression in B cells and DCs was signicantly associated with OS. Furthermore, BTNL9 expression has signicant positive CARE scores. Conclusions: Low BTNL9 expression can prevent the inltration of naïve B cells and DCs in the tumor microenvironment and worsen the outcome in LUAD patients. Besides, these ndings also suggest the potential role of BTNL9 as a prognostic biomarker and a new immuno-target.

BTN or BTNL molecules that meet the signi cant difference between LUAD and adjacent tissue were only BTNL8 and BTNL9. Then survival analysis and select the molecule with signi cant differences, and the result is only BTNL9 in LUAD. BTNL9 was con rmed through the survival analysis of the three databases KM plotter (http://www.kmplot.com) [11], UALCAN (http://ualcan.path.uab.edu/index.html) [12] and OncoLnc (http://www.oncolnc.org/) [13], and subsequent network analysis including immune in ltrates were presented. GEPIA was used to analyze the correlation between BTNL9 and methyltransferase expression, B, and DC cell immune in ltration.

Oncomine database analysis
The mRNA expression of BTNL9 in different tumors was identi ed using the Oncomine database (https://www.oncomine.org/) [14]. The default setting parameters of the database were used, with the threshold parameters set as a P-value of 0.01, a fold change of 1.5, and a top 10% gene ranking.
TIMER database analysis TIMER is a comprehensive resource for systematical analysis of immune in ltrates across diverse cancer types (http://timer.compgenomics.org/) [10,19]. Using Gene_DE module, BTNL9 mRNA expression was calculatedin cross-carcinoma(*: P-value < 0.05; **: P-value <0.01; ***: P-value <0.001). The expression of BTNL9 in LUAD and its correlation with the six immune in ltration cells, namely B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells, used as the Gene and Survival module. Gene_Corr module was used to analyze the correlation between BTNL9 expression and B and DC cells [20].

StarBase database analysis
We used three databases: miRMap (https://mirmap.ezlab.org/) [21], TargetScan (http://www.targetscan.org/vert_72/) [22], and miRWalk (http://mirwalk.umm.uni-heidelberg.de/) [23] to predict the miRNAs that could bind to BTNL9, consider the intersection, verify in the Starbase database (http://starbase.sysu.edu.cn/panCancer.php) [24], and select miRNAs with signi cant differences between co-expression and overall survival. Sangerbox database analysis Sangerbox (http://sangerbox.com/Index) was developed by Hangzhou Mugu Technology Co., Ltd. It mainly uses public data or existing data for immune-related analysis and directly issues a set of charts and analysis reports. Sangerbox uses the R software ESTIMATE package to analyze the immune score and stromal score of each TCGA tumor sample, and GSEA to analyze the enrichment for the KEGG and HALLMARK pathways in the BTNL9 high and low expression groups. Analysis of the relationship between BTNL9 and TMB was determined using Spearman correlation analysis.

Statistical Analysis
The boxplots of BTNL9 expression were extracted from GEPIA. Overall survival was determined from GEPIA, KM plotter, UALCAN, and OncoLnc databases. The expression of BTNL9 across cancer species was extracted in TIMER. RFS and OS HR were extracted from TissGDB database. The correlation between the miRNA and BTNL9 expression was determined in Starbase, as well as the expression of miRNA on the survival of LUAD. BTNL9 expression and single cancer cell biological behavior of LUAD was performed by Pearson correlation analysis. The correlation between BTNL9 expression, the 6 immune cells, and survival analysis was assessed in TIMER. In all the above analyses, P<0.05 was consideredstatistically signi cant.

Results
BTNL9 mRNA expression in LUAD The "Tidyverse" and "Table1" 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 signi cant 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 signi cantly different in LUAD, and both tumor tissues were signi cantly lower than normal tissues ( Figure 1A). BTNL9 signi cantly correlated with the clinical stage, N stage and P53 mutation, while BTNL8 signi cantly correlated with the clinical stage and N stage ( Figure 1B). However, survival analysis revealed that BTNL8 was not signi cantly correlated with OS, while BTNL9 was signi cantly correlated with OS in LUAD ( Figure 1C). PrognoScan [25] analysis of the NSCLC dataset found that BTNL9 expression was signi cantly correlated with RFS and OS in the GSE31210 study, and BTNL9 expression was signi cantly correlated with OS in GSE3141 study (Supplementary Table 1). KM plotter, UALCAN, and OncoLnc were used to validate the ndings, and the results showed that high expression of BTNL9 signi cantly improved the OS of LUAD ( Figure 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 signi cantly 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 signi cantly increased in brain and CNS cancer, four datasets of colorectal cancer, esophageal cancer, leukemia, and lymphoma ( Figure 2A and Supplementary Table 2). We further examined BTNL9 expression using TCGA RNA sequencing data (TIMER). As shown in Figure
Moreover, we analyzed the protein network that binds to BTNL9 through the STRING [29] database and found that there were 7 binding proteins ( Figure 3F). The top 2 binding proteins were predicted using the cytoHubba module of Cytoscape [30] ( Figure 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 (Con dence level is high, and score=0.805), and there is one potential E3 recognizing domain and two potential E3 recognizing motif ( Figure 3H, 3I).

Low BTNL9 expression signi cantly 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 signi cantly enriched in vascular_smooth_muscle_contraction, phosphatidylinositol_signaling_system, and abc_transporters ( Figure 4A). The top 4 KEGGs pathways associated with low expression of BTNL9 were parkinsons_disease, oxidative_phosphorylation, DNA replication, and proteasome ( Figure 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 ( Figure 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 ( Figure 4E

Correlation between BTNL9 and markers of in ltrating 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 signi cantly negatively correlated with TMB (P=1.4E-9) ( Figure 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% ( Figure 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 signi cantly positively correlated with ImmuneScore (r=0.129, P=0.003) and ESTIMATEScore (r=0.106, P=0.016). However, StromalScore was not signi cantly different ( Figure 5C, 5D,5E).
DC also have a trend of difference, and high DC expression of BTNL9 signi cantly prolonged OS in LUAD ( Figure 5F, 5G). We conducted a detailed analysis of TME in ltrated 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 signi cantly 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 signi cantly 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 signi cantly 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.
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 pro le of the cancer sample and the CARE scoring vector, the patient is predicted as a responder or a non-responder [37]. BTNL9 expression has signi cant 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 ( Figure 6, and Table 3).

Discussion
This study systematically analyzed the correlation between human BTNL9 and LUAD and normal adjacent tissues and revealed that BTNL9 can be considered the key gene of butyrophilins in LUAD. Through multiple database survival analysis, this study demonstrated that high BTNL9 expression signi cantly prolonged OS of LUAD. BTNL9 is signi cantly related to B cells, CD4 + T, and macrophages in the TME of LUAD, while high expression of BTNL9 in B and DC cells signi cantly prolongs the OS of LUAD. GEPIA correlation analysis shows that BTNL9 is related to DC cells and naïve B cells but not plasma B cells in LUAD. Moreover, BTNL9 expression is associated with many targeted therapy response.
The low expression of BTNL9 in LUAD is regulated by unknown mechanisms. DNA methylation is the most common form of DNA modi cation, which plays an important role in normal cell physiology, increased DNA methylation, and loss of demethylation, which frequent in cancer. DNMTs play an important role during abnormal DNA methylation. Gene body hypermethylation may activate oncogenes, and promoter hypermethylation often causes the suppression of tumors [38].GEPIA analysis revealed that BTNL9 and DNMTs were signi cantly related to normal lung tissues, but not to LUAD (Fig. 3A, 3B). However, the molecular mechanism of DNMTs regulating BTNL9 in LUAD is unclear. miRNA and lncRNA are common non-coding RNAs, involved in tumor-promoting and suppressing, depending on the type of tumor [38]. We predicted that the 3 miRNAs (hsa-miR-30b-3p, hsa-miR-4709-3p, and hsa-miR-6514-39) were signi cantly positively correlated with BTNL9 in LUAD, and their high expression signi cantly prolongs OS (Fig. 3C, 3D), by exerting anti-cancer activity in LUAD.
Hsu Y-L et al. reported that BTNL9 may act as a tumor suppressor in LUAD and may be regulated by hsa-miR-183-5p, but the speci c regulatory network was not explained [39]. We further predicted that lncRNA AP001462.6 could bind to BTNL9, and the high expression signi cantly prolonged the OS of LUAD (Fig. 3E). Moreover, we analyzed the protein interaction network of BTNL9 and found that the interacting proteins played a role in immune regulation, protease hydrolysis, serine/threonine kinase regulation, etc. (Fig. 3F, 3G), and protease hydrolysis is related to ubiquitination degradation. The analysis also revealed that BTNL9 has a potential E3 recognizing domainbinding site (Fig. 3H, 3I). These indicates that BTNL9 in LUAD may be regulated by DNA methylation and non-coding RNA, and the protein may also be regulated by ubiquitination degradation after translation.
To understand the biological function of BTNL9 in LUAD, GSEA analysis was performed and revealed that low expression of BTNL9 was enriched in energy metabolism (oxidative_phosphorylation, glycolysis, myc_targets v1 [40], and mTORC1_signaling [41]), DNA replication, and protease hydrolysis (Fig. 4A-4D). Metabolic reprogramming triggers selective gene ampli cation and a large gene family which drives cellular functions to promote cancer cell growth and proliferation [40]. The above functions were subsequently veri ed from a single cell perspective and found that BTNL9 expression was signi cantly negatively correlated with cancer cell malignant behaviors such as proliferation, invasion, EMT, metastasis, and hypoxia, etc. in LUAD (Fig. 4E, 4F). This indicates that BTNL9 may play an important and extensive tumor suppressor effect in LUAD.
We further analyzed the relationship between BTNL9 and TME and found that the mutation frequency of BTNL9 in LUAD was about 1.14%, and BTNL9 was signi cantly negatively correlated with TMB (Fig. 5A, 5B). BTNL9 was signi cantly positively correlated with ImmuneScore and ESTIMATEScore (Fig. 5C, 5D, 5E). Previous studies have shown that high ImmuneScore and ESTIMATEScore signi cantly improve the prognosis of LUAD [42]. This shows that BTNL9 plays an important role in TME immune regulation. Moreover, we analyzed the correlation between BTNL9 and TILs and found that BTNL9 was signi cantly negatively correlated with tumor purity, while previous studies have found that low tumor purity is associated with poor prognosis [43]. Although BTNL9 is signi cantly related to B, CD4 + T, and macrophages, survival analysis showed that BTNL9 was only signi cantly related to B and DC cells (Fig. 5F, 5G).
DCs play a pivotal role in shaping the innateand adaptive immune response because they have the unique ability to initiate T cell responses and differentiate them into effector lineages [35], while B cells have antigen presentation, cytotoxicity, and antibody production functions, which are an important part of adaptive immunity [44]. TIMER combined with GEPIA database analysis showed that the expression of BTNL9 on DCs (cDC1s and cDC2s) was not correlated with normal adjacent tissues, but was signi cantly correlated with LUAD tissues (Table 2). cDC1s can migrate to tumor-draining lymph nodes, activate and attract T cells, secrete cytokines, and present antigens in TME, thereby enhancing the function of local cytotoxic T cells [45]. cDC2s present antigens to MHC , activate CD4 + T cells, and effectively polarize TILs into anti-tumor T helper cell 1 (Th1) or Th17 phenotype [46]. BTNL9 on B cells (naïve B cells) was not correlated with normal adjacent tissues but was signi cantly correlated with LUAD tissues, except for plasma B cells, suggesting that BTNL9 regulates the function of naïve B cells in TME. Previous studies have also shown that naïve B cells are down-regulated in advanced NSCLC and are related to poor prognosis [36]. Furthermore, CARE database analysis showed that BTNL9 expression is to be associated with effective target therapy response (Fig. 6, and Table 3).

Conclusions
In conclusion, this study shows that BTNL9 is lowly expressed in LUAD, and the low expression is associated with poor prognosis. BTNL9 may be regulated by both epigenetic regulation and post-transcriptional modi cation. These regulatory mechanisms may change the TME and result in a poor prognostic phenotype and drug resistance. Therefore BTNL9 may play a pivotal role in anti-immune suppression by affecting the TME and cab be considered as a prognostic marker and new immuno-target for LUAD.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Yes.

Availability of data and materials
All data in this study comes from public and open source databases.
Author contributions ZC, TD contributed to the conception and design of the study. WM, JL, ZC, NH, and SZ contributed to data collection. ZC, WM, and TD analyzed and interpreted the data. ZC drafted the manuscript. All authors approved the nal version of the report.

Con icts of interest
The authors declare no potential con ict of interest.