Skip to main content
Download PDF

Comprehensive analysis of BTNL9 as a prognostic biomarker correlated with immune infiltrations in thyroid cancer

BMC Medical Genomics volume 16, Article number: 234 (2023) Cite this article

Abstract

Background

Thyroid cancer (THCA) is the most common type of endocrine cancers, and the disease recurrences were usually associated with the risks of metastasis and fatality. Butyrophilin-like protein 9 (BTNL9) is a member of the immunoglobulin families. This study investigated the prognostic role of BTNL9 in THCA.

Methods

Gene enhancers of BTNL9 were identified by interrogating H3K27ac ChIP-seq and RNA-seq data of papillary thyroid cancer (PTC) and benign thyroid nodule (BTN) tissues. Meanwhile, BTNL9 expression level was verified by qRT-PCR in 30 pairs of primary THCA and adjacent normal tissues. Clinicopathological and RNA sequencing data were obtained from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) to analyze the relations between BTNL9 expression and immune cell infiltration, chemokines/cytokines, immune checkpoint genes, clinical parameters and prognosis values. Besides, survival analysis combining BTNL9 expression and immune cell infiltration scores was conducted. Functional enrichment analysis was performed to investigate the potential biological mechanisms. Cox regression analyses were used to explore independent clinical indicators, and a nomogram model incorporating BTNL9 expression with clinical parameters was established.

Results

BTNL9 showed significantly stronger H3K27ac modifications in BTN than PTC tissues at the promoter region (chr5: 181,035,673–181,047,436) and gene body (chr5: 181,051,544–181,054,849). The expression levels of BTNL9 were significantly down-regulated in THCA samples compared to normal tissues, and were strongly associated with different tumor stages, immune cell infiltrations, chemokines/cytokines and immune checkpoint genes in THCA. Functional enrichment analyses indicated that BTNL9 was involved in immune-related and cancer-related pathways. The Kaplan–Meier analysis showed lower BTNL9 expression was associated with poorer progression-free interval (PFI). BTNL9 expression and pathologic stages were independent prognostic indicators of PFI in THCA.

Conclusions

The results implied an important role of BTNL9 in the tumor progression, with the possibility of serving as a novel prognostic biomarker and a potential therapeutic target for THCA.

Peer Review reports

Background

Thyroid cancer (THCA) is the most common type of endocrine cancers, taking up about 90% of all endocrine tumors and 2% of all systemic malignancies [1, 2], and it ranked the fifth most prevalent tumor in women [3]. Generally, THCA could be categorized into different subtypes based on histologic features, among which papillary thyroid cancer (PTC) accounts for over 80% of the cases [2]. About one third of THCA patients might evolve with disease recurrences, facing the risks of metastasis and THCA-related fatality [4, 5]. Moreover, with the development of imaging technology and increased use of neck ultrasound for the routine detection of thyroid neoplasms, the thyroid nodules were detected in approximately 50% adults [6]. Thus, it is of great significance to differentiate between malignant and benign nodules so as to reduce overdiagnosis and unnecessary treatment for people with thyroid nodules. However, the fine-needle aspiration guided by ultrasound could only help to make differentiation diagnosis in about 75% nodules [7, 8], and the remaining indeterminate nodules hindered the medical management of these patients. Therefore, it is of clinical significance to find new molecular biomarkers for the differentiation and prognosis evaluation.

Butyrophilin-like protein 9 (BTNL9) is a member of the immunoglobulin butyrophilin and butyrophilin-like (BTNL) families, modulating the T cell response and impacting inflammatory disorders and cancers [9]. It is reported that recombinant BTNL9–Fc could bind with activated T cells, B cells, dendritic cells and macrophages [10, 11]. In immune cells, BTNL9 was mostly expressed in B cells yet its functions were not well understood. Few studies reported the roles and mechanisms of BTNL9 in cancer. BTNL9 was reported to have lower expression levels in various cancers compared to normal tissues, including colon cancer, lung adenocarcinoma, osteosarcoma, breast cancer and uveal melanoma [12,13,14,15,16]. However, the prognosis value and immunology of BTNL9 in THCA have not been elucidated.

In this study, the H3K27ac modification and RNA expression differences of BTNL9 gene between PTC and benign thyroid nodules (BTN) were characterized, revealing the epigenetic modifications on gene transcription spectrum. The expression level of BTNL9 was further validated in 30 pairs of THCA and corresponding normal tissues. Furthermore, comprehensive analyses were performed including the correlation between BTNL9 expression and clinical parameters, prognostic value, immune cell infiltration, immune biomarkers and immune checkpoints, with functional enrichment analyzed on differentially expressed genes (DEGs) to explore the functions and potential mechanisms of BTNL9, raising the possibility of BTNL9 serving as a novel prognostic indicator and an immunotherapeutic target in THCA.

Methods

Data collection and analysis

Gene mRNA expression data with clinical information was downloaded from the TCGA database (https://portal.gdc.cancer.gov/) and the GTEx projects (Genotype-Tissue Expression, https://gtexportal.org/). Statistical analysis was performed using R software (v. 3.6.3). The differential expression analysis of BTNL9 gene was analyzed with GEPIA online database (Gene Expression Profiling Interactive Analysis) (http://gepia.cancer-pku.cn/), which is a comprehensive platform with TCGA and GTEx data analyzed for cancer research [17].

Survival prognosis analysis

The Kaplan-Meier plotter (http://kmplot.com), an online interactive platform, which could estimate the effects of gene expression on prognosis in various cancers [18], was used to investigate the correlation between BTNL9 expression and prognosis in THCA. Besides, subgroup survival analysis was performed using different clinical parameters. The hazard ratio (HR) and log-rank p-value were estimated.

Cox proportional hazards models were performed to evaluate prognostic indicators of THCA. Univariate and multivariate Cox regression analyses were performed with survival package (v. 3.2–10). Nomogram model incorporating BTNL9 expression and clinical indicators were generated using rms (v. 6.2-0) and survival package.

Immune infiltration analysis

The Tumor Immune Estimation Resource (TIMER2.0) (http://timer.cistrome.org/) [19], a comprehensive online platform for tumor-infiltrating assessment, was used to analyze the immune cell infiltrations in THCA. The associations between BTNL9 gene and various chemokines/cytokines were evaluated via TIMER database.

The relationships between BTNL9 expression and infiltration levels of 24 types of immune cells [20] were analyzed via ssGSEA algorithm using GSVA package (v. 1.34.0) [21]. The associations between BTNL9 expression and immune checkpoint genes were visualized by heatmap plots using ggplot2 R package (v. 3.3.3). Kaplan-Meier survival analysis with the combined data of BTNL9 expression and immune cell infiltration scores was conducted with survminer (v. 0.4.9) and survival (v. 3.2–10) packages.

Gene enrichment analysis

The Linked Omics database (http://www.linkes.org/login.php) is a web-based tool for the analysis of multi-omics data across different cancer types [22]. The genes associated with BTNL9 were analyzed and visualized by the Linked Omics.

The DEGs were analyzed by DESeq2 package (v.1.26.0) with the median BTNL9 expression level set as cut-off value to define low and high expression subgroups. DEGs with |Log2(Fold Change)| >1 and adjusted p<0.05 were then subjected to functional enrichment analysis. The Gene Ontology (GO) terms annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed with ClusterProfiler package (v. 3.14.3) [23]. Gene set enrichment analysis (GSEA) [24] using predefined gene sets from MSigDB Collections (https://www.gsea-msigdb.org/gsea/msigdb/collections.jsp) was conducted, and results with false discovery rate (FDR) < 0.25 and adjusted p<0.05 were defined as significant. The SRTING Database (https://cn.string-db.org/), an online tool for protein-protein interaction networks analysis, was used to investigate the integrated network of BTNL9.

H3K27ac ChIP-seq and RNA-seq data

The detailed protocols of ChIP-seq and RNA-seq were reported in our previous article [25]. The data can be found at Genome Sequence Archive in National Genomics Data Center, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (Accession number HRA000779).

Real-time quantitative PCR analysis

30 pairs of excised THCA and paracancerous thyroid tissues were obtained from the First Affiliated Hospital of Sun Yat-sen University after institutional review board approval and informed patient consent. Total mRNAs of tissues were extracted using RNAsimple Total RNA Kit (#DP419; TIANGEN Biotech). 500ng mRNA was used as template for reverse transcription with PrimeScript RT Master Mix (RR036A; Takara), and SYBR Green Master Mix (A25742; Thermo Fisher Scientific) was used to detect cDNA amplification. GAPDH was used to normalize gene expression. The primers for BTNL9 were as follows: forward 5’- ATGGTGGACCTCTCAGTCTCC-3’, reverse 5’-GCCAGGATGGGATACTCAGG-3’. GAPDH primers: forward 5’- ACTTCAACAGCGACACCCACTC-3’, reverse 5’- TCTCTTCCTCTTGTGCTCTTGCT-3’. Relative expression of genes was calculated using the 2− ΔΔCT method.

Statistics

Statistical analysis was performed using R software (v 3.6.3). One-way ANOVA test and student t test were adopted in comparing the clinical features of high-expression and low-expression groups; Chi-Square test, Wilcoxon and Kruskal-Wallis tests were used when needed. Spearman’s correlation coefficient was used to evaluate the correlations with immune cell infiltrations. Statistical significance was defined as p < 0.05.

Results

Differentiated BTNL9 expression between thyroid cancer and normal tissues

A flow chart illustrating the analysis procedure of the study is shown in Fig. 1. H3K27ac ChIP-seq and RNA-seq data of PTC and BTN were integrated to map the enhancers of BTNL9 gene. As shown in Fig. 2A, BTN samples had stronger H3K27ac signals than PTC samples at the promoter (chr5: 181,035,673–181,047,436) and gene body (chr5: 181,051,544–181,054,849) of BTNL9. The RNA-seq data revealed that BTNL9 expression was higher in BTN than PTC samples. Additionally, qRT-PCR on 30 pairs of primary samples verified that BTNL9 expression levels were significantly down-regulated in THCA than matched paracancerous tissues frequently (Fig. 2F). Consistently, when analyzing the data from TCGA and GTEx, significantly lower expression levels of BTNL9 in THCA than normal controls (Fig. 2B) and matched normal samples (n = 58) were observed (Fig. 2D). The differences in epigenetic H3K27ac modification and expression indicated that BTNL9 might be a potential diagnostic biomarker for the differentiation between THCA and benign nodules.

Fig. 1
figure 1

Flow chart shows the analysis procedure of the study. PTC, papillary thyroid carcinoma; BTN, benign thyroid nodule; TCGA, The Cancer Genome Atlas; GTEx, Genotype-Tissue Expression; THCA, thyroid carcinoma

Full size image
Fig. 2
figure 2

Expression levels of BTNL9 gene in THCA. (A) Gene tracks of BTN-specific enhancers (chr5:181,035,673–181,047,436; chr5:181,051,544–181,054,849). The adjacent gene BTNL9 showed much stronger signals in H3K27ac ChIP-seq and RNA-seq tracks in BTN tissues than PTC tissues. (B) The expression level of BTNL9 was significantly lower in THCA compared with normal samples based on TCGA and GTEx datasets. (C)BTNL9 expression in different cancer types from TCGA and GTEx datasets. (D) Decreased expression of BTNL9 in THCA compared with the matched normal tissues from TCGA datasets. (E) Higher level of BTNL9 transcription is associated with more favorable PFI in patients with THCA. Analysis was performed with TCGA (THCA) datasets. (F) Relative BTNL9 mRNA expression levels in 30 pairs of THCA and paired paracancerous tissues. Expression analysis was normalized against GAPDH. *p < 0.05, **p < 0.01, ***p < 0.001. THCA, thyroid carcinoma; PTC, papillary thyroid carcinoma; BTN, benign thyroid nodule; PFI, progress-free interval; TCGA, The Cancer Genome Atlas; GTEx, Genotype-Tissue Expression. ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical and endocervical cancers; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; DLBC, lymphoid neoplasm diffuse large B-cell lymphoma; ESCA, esophageal carcinoma; GBM, glioblastoma multiforme; HNSC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; KIRC, kidney renal clear cell carcinoma; KIRP, kidney renal papillary cell carcinoma; LAML, acute myeloid leukemia; LGG, brain lower grade glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PRAD, prostate adenocarcinoma; READ, rectum adenocarcinoma; SARC, sarcoma; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; TGCT, testicular germ cell tumors; THYM, thymoma; UCEC, uterine corpus endometrial carcinoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma

Full size image

BTNL9 expression levels in other 32 types of cancer were analyzed with TCGA and GTEx datasets (Fig. 2C). It revealed that BTNL9 expression was significantly down-regulated in various cancers, including adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical and endocervical cancer (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), testicular germ cell tumor (TGCT), uterine corpus endometrial carcinoma (UCEC), and uterine carcinosarcoma (UCS).

Prognostic analysis and clinical relevance of BTNL9 expression in thyroid cancer

Kaplan–Meier survival curves demonstrated that high BTNL9 expression was correlated with favorable progression free interval (PFI) (Fig. 2E) and relapse-free survival (RFS) (Fig. 3A) in patients with THCA based on TCGA datasets. Specifically, subgroup survival analysis showed that increased BTNL9 expression was associated with favorable RFS in female patients (Fig. 3B), tumor with low neoantigen load (Fig. 3C) and patients with stage 1 THCA (Fig. 3D). In stage 4 THCA subgroup, greater BTNL9 expression was related to better overall survival (OS) outcomes (Fig. 3E). These results suggested the potential of BTNL9 serving as a prognostic biomarker in THCA.

Fig. 3
figure 3

Kaplan-Meier survival curves for patients with THCA from KM plotter web-based tool. (A) Higher expression of BTNL9 was associated with better RFS in patients with THCA. (B) Subgroup analysis indicated that BTNL9 expression was associated with more favorable RFS in female THCA patients. (C)BTNL9 expression was associated with more favorable RFS in THCA patients with low neoantigen load. (D) Higher expression of BTNL9 was associated with better RFS in patients with Stage 1 THCA. (E) Higher expression of BTNL9 was associated with better OS in patients with Stage 4 THCA.

RFS, relapse free survival; OS, overall survival; THCA, thyroid carcinoma

Full size image

To analyze the relationships between BTNL9 expression and clinical characteristics in patients with THCA, 502 patients from TCGA datasets were categorized into high and low expression groups with the mean value of BTNL9 expression set as dividing threshold (Table 1). It showed that BTNL9 expression was significantly associated with different clinical T stages, N stages, extrathyroidal extension, pathologic stages, histological types, primary neoplasm focus types, residual tumors, thyroid gland disorder history and disease progression (Fig. 4) in THCA patients.

Table 1 Correlation between BTNL9 expression and clinical characteristics in patients with thyroid cancer
Full size table
Fig. 4
figure 4

Association between BTNL9 expression and clinical features in patients with THCA based on TCGA datasets. BTNL9 expression levels were stratified by (A) tumor T stage, (B) tumor N stage, (C) extrathyroidal extension, (D) pathologic stage, (E) histological type, (F) primary neoplasm focus type, (G) residual tumor, (H) thyroid gland disorder history and (I) PFI event. TCGA, The Cancer Genome Atlas; THCA, thyroid carcinoma; PFI, progress-free interval. *p < 0.05, **p < 0.01, ***p < 0.001

Full size image

Correlation between BTNL9 expression and immune cell infiltration

The tumor microenvironment components played a vital role in THCA, exerting both antitumor and tumor-promoting functions [26,27,28,29,30,31]. The enrichment levels of 24 types of immune cells were assessed between groups of high and low BTNL9 expression, and it revealed that T cells, Treg cells (regulatory T cell), Th1 cells (T helper 1 cell), Th2 cells (T helper 2 cell), T helper cells, Tem cells (effector memory T cell), Tcm cells (central memory T cell), neutrophils, mast cells, macrophages, eosinophils, DC (dendritic cells), iDC (interdigitating cell), aDC (activated dendritic cells), cytotoxic cells and B cells were significantly enriched in low BTNL9 expression group (Fig. 5A).

Fig. 5
figure 5

Association between BTNL9 expression and immune infiltrations in THCA. (A) Box plot illustrating the enrichment of 24 subtypes of immune cells in high and low BTNL9 expression groups. (B) Lollipop chart showing the correlation between BTNL9 expression and different immune cells. (C)BTNL9 copy number variations were related to the infiltration levels of B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils and dendritic cells in THCA. (D) Correlations between BTNL9 expression and infiltration levels of different immune cells. (E) Heatmap indicating correlations between BTNL9 expression and immune checkpoint genes. THCA, thyroid carcinoma; Tcm cell, central memory T cell; Treg cell, regulatory T cell; Th1 cell, T helper 1 cell; Th2 cell, T helper 2 cell; Tem cell, effector memory T cell; TFH cell, follicular helper T cell; NK cell, natural killer cell; DC, dendritic cells; pDC, plasmacytoid dendritic cells; iDC, interdigitating cell; aDC, activated dendritic cells. *p < 0.05, **p < 0.01, ***p < 0.001

Full size image

Besides, it revealed that BTNL9 expression was positively correlated with pDC (plasmacytoid dendritic cells, r = 0.343, p < 0.001), Tgd cells (r = 0.098, p = 0.028) and NK cells (natural killer cells, r = 0.095, p = 0.032). In contrast, negative correlations were found between BTNL9 expression and enrichment levels of NK CD56dim cells (r =-0.092, p = 0.038), Tcm cells (r = -0.101, p = 0.024), TFH cells (follicular helper T cells, r = -0.129, p = 0.004), eosinophils (r = -0.153, p < 0.001), Tem cells (effector memory T cells, r = -0.155, p < 0.001), T helper cells (r = -0.21, p < 0.001), cytotoxic cells (r = -0.282, p < 0.001), mast cells (r = -0.305, p < 0.001), B cells (r = -0.342, p < 0.001), T cells (r = -0.378, p < 0.001), neutrophils (r = -0.485, p < 0.001), Th2 cells (r = -0.529, p < 0.001), iDC (r = -0.53, p < 0.001), Th1 cells (r = -0.57, p < 0.001), macrophages (r = -0.598, p < 0.001), aDC (r = -0.597, p < 0.001), DC (r = -0.622, p < 0.001), Treg cells (r = -0.641, p < 0.001) (Fig. 5B,D).

BTNL9 expression showed strong relationships with a large proportion of chemokines/cytokines (Additional file 1: Table S1), which indicated that BTNL9 expression might get involved in the immune-related biological processes in coordination with these genes. Copy number variations of BTNL9 gene were found to have associations with the infiltration levels of B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils and dendritic cells in THCA (Fig. 5C). The associations between BTNL9 expression and common immune checkpoint genes were explored in THCA, and most of the immune checkpoint genes were significantly negatively correlated with BTNL9 expression (Fig. 5E).

Furthermore, the data of BTNL9 expression combined with immune cell infiltration from TCGA dataset was explored by Kaplan-Meier analysis. It demonstrated that in THCA, patients with high BTNL9 expression and low infiltration scores of Tcm cells, Tem cells or pDC showed favorable outcomes of PFI (Fig. 6A-C), whereas those with high BTNL9 expression combined with high scores of Th1 cells, Th2 cells, neutrophils, T cells, Treg cells, iDC, aDC, mast cells, cytotoxic cells or macrophages (Fig. 6D-M) presented better PFI outcomes in THCA patients.

Fig. 6
figure 6

Kaplan-Meier survival curves of PFI stratified by immune cell scores and BTNL9 expression in THCA. PFI curves analyzed with the combination of BTNL9 expression level and (A) Tcm cell scores, (B) Tem cell scores, (C) pDC scores, (D) Th1 cell scores, (E) Th2 cell scores, (F) neutrophil cell scores, (G) T cell scores, (H) Treg cell scores, (I) iDC scores, (J) aDC scores, (K) mast cell scores, (L) cytotoxic cell scores and (M) macrophage cell scores. PFI, progress-free interval; THCA, thyroid carcinoma; Tcm cell, central memory T cell; Treg cell, regulatory T cell; Th1 cell, T helper 1 cell; Th2 cell, T helper 2 cell; Tem cell, effector memory T cell; pDC, plasmacytoid dendritic cells; iDC, interdigitating cell; aDC, activated dendritic cells

Full size image

Functional enrichment analysis of DEGs associated with BTNL9 in THCA

A total of 1595 DEGs associated with BTNL9 were identified in THCA, including 733 up-regulated and 862 down-regulated genes (Additional file 2: Table S2). Besides, the co-expression network of BTNL9 gene was explored by the Linkedomics platform with top 50 positively-related genes and negatively-related genes visualized in the forms of heatmaps (Fig. 7A, B).

Fig. 7
figure 7

Functional enrichment analysis of differentially expressed genes (DEGs) associated with BTNL9 gene in THCA. (A) Top 50 genes positively related with BTNL9 in THCA. (B) Top 50 genes negatively related with BTNL9 in THCA. (C) GO enrichment analysis of co-expressed genes. (D) KEGG pathway analysis [32] of co-expressed genes. (E) BTNL9 interaction network analysis using STRING tool. (F-K) GSEA plots showing DEGs enriched in immune-related pathways and cancer-related pathways. THCA, thyroid carcinoma; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis

Full size image

Enrichment analyses of the DEGs were performed. GO annotations revealed that related genes were mainly involved in immune-related pathways such as leukocyte migration, T cell activation and cell chemotaxis (Fig. 7C). Similarly, KEGG enrichment analysis [32] revealed that these genes were mostly enriched in immune functions such as chemokine signaling pathway, complement and coagulation cascades, and IL-17 signaling pathway (Fig. 7D). The detailed pathway analysis was presented in Additional file 2: Table S2. GSEA analysis showed these genes were enriched in tumorigenesis and immune-related pathways, such as TNFα -signaling via NF-kappa B, P53 pathway, and inflammatory response (Fig. 7F-K). In addition, STRING online platform was used to investigate the protein–protein interaction (PPI) network, and 6 proteins were found to have direct interactions with BTNL9 (Fig. 7E).

Prognostic analysis of clinical indicators in thyroid cancer

Cox regression analysis was performed to study the indicators of PFI in THCA. It showed that in univariate model, pathologic stage, tumor extrathyroidal extension and BTNL9 expression were related to PFI (Fig. 8A) (Table 2). Low BTNL9 expression is associated with higher hazard ratio of PFI in THCA. The factors with statistical significance in univariate analysis were then enrolled into multivariate Cox regression, and it turned out that BTNL9 expression and pathologic stage were independent prognostic indicators (Table 2). The risk score list of the Cox regression analysis was provided in Additional file 3: Table S3. Furthermore, nomogram model was constructed to predict 10-year PFI in THCA patients using the above clinical features, with the calibration curve evaluating the consistency between predicted probability and ideal model plotted (Fig. 8B-C).

Fig. 8
figure 8

Prognostic prediction analysis of THCA. (A) Forest plot for univariate analysis of PFI in THCA based on TCGA datasets. (B) Nomogram plot predicting the probability of 10-year PFI in THCA patients. (C) Calibration plot for 10-year PFI prediction. THCA, thyroid carcinoma; PFI, progress-free interval

Full size image
Table 2 Cox regression analysis of prognostic factors for progression free interval in thyroid cancer
Full size table

Discussion

Cancer biomarkers, which included genetic alterations, epigenetic changes, proteins and pathological features, had been used in facilitating accurate diagnosis and prognostic evaluations in THCA, providing information for precise clinical management. It is recommended by American Thyroid Association guidelines that molecular testing is useful in making appropriate treatment regimen and reducing unnecessary surgical operation [4, 33]. With the development of high-throughput screening, more novel biomarkers were identified. In this study, by integrating H3K27ac ChIP-seq and RNA-seq data, BTNL9 was found to have stronger H3K27ac signals with higher expression in BTN samples than PTC samples. The expression differences were further validated in PTC and corresponding adjacent normal tissues. Besides, high BTNL9 expression in THCA was found to be related with favorable PFI in THCA based on TCGA datasets, and comprehensive bioinformatics analysis was conducted, suggesting the potential value of BTNL9 in differential diagnosis of thyroid nodules and clinical prognosis.

Epigenetic alterations greatly impact cancer development. H3K27ac is an established mark of transcriptional activation. Our data presented epigenetic regulations of BTNL9 gene with differential expression spectrums in PTC and BTN samples. Besides from THCA, decreased expression levels of BTNL9 gene were also presented in various cancers according to TCGA datasets, which is consistent with previous studies [12,13,14,15,16]. Kaplan-Meier analysis suggested the prognostic value of BTNL9 in THCA. Analyses of BTNL9 expression in different clinical and pathologic stages showed decreased expression in more advanced stages, indicating BTNL9 might participate in tumor progression. As for clinical parameters, lower BTNL9 expression was shown in tall cell histological type, unifocal neoplasms, R1 and R2 residual tumors, and lymphocytic thyroiditis compared with other subtypes. The results indicated the potential involvement of BTNL9 in tumorigenesis and progression as the above clinical factors were regarded as established risk factors for THCA [34, 35].

The tumor microenvironment components played a vital role in tumor initiation, progression and metastasis [36]. Numerous immune cells infiltrated the THCA microenvironment, exerting both antitumor and promoting functions, and in some way influence patients’ outcomes [26,27,28,29,30,31]. A number of subsets of immune cells were found to have effects on THCA and be correlated with disease outcomes. It was reported that immature DC subsets functioned as immune suppressor in tumor progression while matured DC showed anti-tumor effects [37]. Treg cells shut down antitumor immune response and promote tumor growth, and increased Treg cells in metastatic lymph nodes were associated with aggressive THCA [38], while in PTC a higher Treg cell infiltration was related to disease progression, extrathyroidal extension and recurrence [38,39,40]. Mast cell infiltration was associated with extrathyroidal extension in PTC [28]. Tumor-associated neutrophil density was correlated with tumor size in THCA samples, indicating the potential of promoting cancer growth [41]. High density of tumor-associated macrophages was related to more malignant biological behaviors, including lymph node metastasis in PTC [42], decreased survival in advanced THCA [27], and distant lung metastasis [31]. BTNL9 was shown to have negative correlations with these immune cells, indicating it might get involved in immunity microenvironment modulation. Immune-related cytokines and chemokines were investigated in THCA, and some of them exerted tumor-promotion effect. CCL2 was reported to promote lymph node metastasis and recurrence in PTC [43], and IL-10 was related to tumor size and invasion in THCA [44]. FOXP3 was found to have associations with DTC aggressiveness and tumor diameter [45]. The correlations between BTNL9 and these cytokines/chemokines indicated its involvement in tumor immunity process. The immune checkpoint PD-1/PD-L1 axis (CTLA-4, TIGIT, TIM3) activation played an important part in thyroid carcinogenesis [46]. PD-1 expression was related to DTC aggressiveness and unfavorable survival outcomes in THCA patients [47, 48]. Previous study showed that PD-1 inhibitor could be used in the treatment of ATC [49], and several ongoing clinical trials are evaluating the effects of combining PD-1 inhibitor with other therapies in advanced THCA treatment. BTNL9 expression was negatively related with many immune checkpoint genes including PD-1 axis, suggesting the possibility of BTNL9 as a potential therapeutic target in THCA.

Furthermore, the functional enrichment analysis showed BTNL9 was involved in the immune-related and tumor-related regulating signaling pathways, indicating that low expression of BTNL9 might promote cancer development in THCA. Cox regression analysis suggested that BTNL9 expression was an independent prognostic indicator in THCA, and a nomogram model incorporating BTNL9 expression was established to predict 10-year PFI, which could provide a clinical reference for identifying high-risk patients. However, some limitations should be pointed out. First, the analyses were mainly performed by bioinformatics method based on public database, and further validation experiments need to be conducted. Second, the mechanisms by which BTNL9 take part in immune regulation remained unclear and warrant further in vivo and in vitro studies.

Conclusion

In summary, the epigenetic H3K27ac modifications and expression levels of BTNL9 gene were decreased in THCA. BTNL9 expression was an independent prognostic indicator of PFI, and low expression was related to unfavorable prognosis in THCA patients. The BTNL9 expression was closely associated with various immune cell infiltrations and immune checkpoints. The DEGs associated with BTNL9 were mainly involved in immune-related and cancer-related pathways. The analyses raised the possibility of BTNL9 being a novel prognostic biomarker and a potential therapeutic target in THCA.

Data Availability

The datasets presented in this study can be found in online repositories. The names of the repositories and accession number can be found in the article. TCGA : https://portal.gdc.cancer.gov/, GTEx: https://gtexportal.org/, GEPIA2: http://gepia.cancer-pku.cn/, Kaplan-Meier plotter: http://kmplot.com,

TIMER2.0: http://timer.cistrome.org/, Linked Omics: http://www.linkes.org/login.php, STRING: https://cn.string-db.org/. H3K27ac ChIP-seq and RNA-seq data can be found in Genome Sequence Archive (GSA): https://ngdc.cncb.ac.cn/gsa/ (Accession number HRA000779).

References

  1. Lim H, Devesa SS, Sosa JA, Check D, Kitahara CM. Trends in thyroid Cancer incidence and mortality in the United States, 1974–2013. JAMA. 2017;317(13):1338–48.

    PubMed  PubMed Central  Google Scholar 

  2. Fagin JA, Wells SJ. Biologic and clinical perspectives on thyroid Cancer. N Engl J Med. 2016;375(11):1054–67.

    CAS  PubMed  PubMed Central  Google Scholar 

  3. Jemal A, Siegel R, Xu J, Ward E. Cancer statistics, 2010. CA Cancer J Clin. 2010;60(5):277–300.

    PubMed  Google Scholar 

  4. Haugen BR, Alexander EK, Bible KC, Doherty GM, Mandel SJ, Nikiforov YE, Pacini F, Randolph GW, Sawka AM, Schlumberger M, et al. 2015 american thyroid Association Management Guidelines for adult patients with thyroid nodules and differentiated thyroid Cancer: the american thyroid Association Guidelines Task Force on thyroid nodules and differentiated thyroid Cancer. THYROID. 2016;26(1):1–133.

    PubMed  PubMed Central  Google Scholar 

  5. Gharib H, Papini E, Garber JR, Duick DS, Harrell RM, Hegedüs L, Paschke R, Valcavi R, Vitti P, AMERICAN ASSOCIATION OF CLINICAL ENDOCRINOLOGISTS, AMERICAN COLLEGE OF ENDOCRINOLOGY, AND ASSOCIAZIONE MEDICI ENDOCRINOLOGI MEDICAL GUIDELINES FOR CLINICAL PRACTICE FOR THE DIAGNOSIS AND MANAGEMENT OF THYROID NODULES. –2016 UPDATE. ENDOCR PRACT. 2016;22(5):622–39.

    PubMed  Google Scholar 

  6. Jiang H, Tian Y, Yan W, Kong Y, Wang H, Wang A, Dou J, Liang P, Mu Y. The prevalence of thyroid nodules and an analysis of related lifestyle factors in Beijing Communities. Int J Environ Res Public Health. 2016;13(4):442.

    PubMed  PubMed Central  Google Scholar 

  7. Baloch ZW, LiVolsi VA, Asa SL, Rosai J, Merino MJ, Randolph G, Vielh P, DeMay RM, Sidawy MK, Frable WJ. Diagnostic terminology and morphologic criteria for cytologic diagnosis of thyroid lesions: a synopsis of the National Cancer Institute thyroid fine-needle aspiration state of the Science Conference. DIAGN CYTOPATHOL. 2008;36(6):425–37.

    PubMed  Google Scholar 

  8. Bongiovanni M, Spitale A, Faquin WC, Mazzucchelli L, Baloch ZW. The Bethesda System for reporting thyroid cytopathology: a meta-analysis. Acta Cytol. 2012;56(4):333–9.

    PubMed  Google Scholar 

  9. Melandri D, Zlatareva I, Chaleil R, Dart RJ, Chancellor A, Nussbaumer O, Polyakova O, Roberts NA, Wesch D, Kabelitz D, et al. The γδTCR combines innate immunity with adaptive immunity by utilizing spatially distinct regions for agonist selection and antigen responsiveness. NAT IMMUNOL. 2018;19(12):1352–65.

    CAS  PubMed  PubMed Central  Google Scholar 

  10. Arnett HA, Escobar SS, Viney JL. Regulation of costimulation in the era of butyrophilins. Cytokine. 2009;46(3):370–5.

    CAS  PubMed  Google Scholar 

  11. Yamazaki T, Goya I, Graf D, Craig S, Martin-Orozco N, Dong C. A butyrophilin family member critically inhibits T cell activation. J IMMUNOL. 2010;185(10):5907–14.

    CAS  PubMed  Google Scholar 

  12. Ho XD, Phung P, Q LV, Reimann HNV, Prans E, Kõks E, Maasalu G, Le K. Whole transcriptome analysis identifies differentially regulated networks between osteosarcoma and normal bone samples. Exp Biol Med (Maywood). 2017;242(18):1802–11.

    CAS  PubMed  Google Scholar 

  13. Lebrero-Fernández C, Wenzel UA, Akeus P, Wang Y, Strid H, Simrén M, Gustavsson B, Börjesson LG, Cardell SL, Öhman L, et al. Altered expression of Butyrophilin (BTN) and BTN-like (BTNL) genes in intestinal inflammation and colon cancer. Immun Inflamm Dis. 2016;4(2):191–200.

    PubMed  PubMed Central  Google Scholar 

  14. Jiang Z, Liu F. Butyrophilin-Like 9 (BTNL9) suppresses Invasion and correlates with favorable prognosis of Uveal Melanoma. Med Sci Monit. 2019;25:3190–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  15. Bao Y, Wang L, Shi L, Yun F, Liu X, Chen Y, Chen C, Ren Y, Jia Y. Transcriptome profiling revealed multiple genes and ECM-receptor interaction pathways that may be associated with breast cancer. CELL MOL BIOL LETT. 2019;24:38.

    PubMed  PubMed Central  Google Scholar 

  16. Fang J, Chen F, Liu D, Gu F, Chen Z, Wang Y. Prognostic value of immune checkpoint molecules in breast cancer. Biosci Rep 2020, 40(7).

  17. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  18. Nagy Á, Munkácsy G, Győrffy B. Pancancer survival analysis of cancer hallmark genes. Sci Rep. 2021;11(1):6047.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: a web server for Comprehensive Analysis of Tumor-Infiltrating Immune cells. CANCER RES. 2017;77(21):e108–10.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, Angell H, Fredriksen T, Lafontaine L, Berger A, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013;39(4):782–95.

    CAS  PubMed  Google Scholar 

  21. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Vasaikar SV, Straub P, Wang J, Zhang B. LinkedOmics: analyzing multi-omics data within and across 32 cancer types. NUCLEIC ACIDS RES. 2018;46(D1):D956–63.

    CAS  PubMed  Google Scholar 

  23. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.

    CAS  PubMed  PubMed Central  Google Scholar 

  25. Zhang L, Xiong D, Liu Q, Luo Y, Tian Y, Xiao X, Sang Y, Liu Y, Hong S, Yu S, et al. Genome-wide histone H3K27 acetylation profiling identified genes correlated with prognosis in papillary thyroid carcinoma. Front Cell Dev Biol. 2021;9:682561.

    PubMed  PubMed Central  Google Scholar 

  26. Liotti F, Visciano C, Melillo RM. Inflammation in thyroid oncogenesis. AM J CANCER RES. 2012;2(3):286–97.

    PubMed  PubMed Central  Google Scholar 

  27. Ryder M, Ghossein RA, Ricarte-Filho JC, Knauf JA, Fagin JA. Increased density of tumor-associated macrophages is associated with decreased survival in advanced thyroid cancer. Endocr Relat Cancer. 2008;15(4):1069–74.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. Melillo RM, Guarino V, Avilla E, Galdiero MR, Liotti F, Prevete N, Rossi FW, Basolo F, Ugolini C, de Paulis A, et al. Mast cells have a protumorigenic role in human thyroid cancer. Oncogene. 2010;29(47):6203–15.

    CAS  PubMed  Google Scholar 

  29. Cho JW, Kim WW, Lee YM, Jeon MJ, Kim WG, Song DE, Park Y, Chung KW, Hong SJ, Sung TY. Impact of tumor-associated macrophages and BRAF(V600E) mutation on clinical outcomes in patients with various thyroid cancers. Head Neck. 2019;41(3):686–91.

    PubMed  Google Scholar 

  30. Jung KY, Cho SW, Kim YA, Kim D, Oh BC, Park DJ, Park YJ. Cancers with higher density of Tumor-Associated Macrophages were Associated with Poor Survival Rates. J Pathol Transl Med. 2015;49(4):318–24.

    PubMed  PubMed Central  Google Scholar 

  31. Li XJ, Gangadaran P, Kalimuthu S, Oh JM, Zhu L, Jeong SY, Lee SW, Lee J, Ahn BC. Role of pulmonary macrophages in initiation of lung metastasis in anaplastic thyroid cancer. INT J CANCER. 2016;139(11):2583–92.

    CAS  PubMed  Google Scholar 

  32. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. NUCLEIC ACIDS RES. 2000;28(1):27–30.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. Patel KN, Yip L, Lubitz CC, Grubbs EG, Miller BS, Shen W, Angelos P, Chen H, Doherty GM, Fahey TR, et al. The American Association of Endocrine Surgeons Guidelines for the definitive Surgical management of thyroid disease in adults. ANN SURG. 2020;271(3):e21–e93.

    PubMed  Google Scholar 

  34. Kondo T, Ezzat S, Asa SL. Pathogenetic mechanisms in thyroid follicular-cell neoplasia. NAT REV CANCER. 2006;6(4):292–306.

    CAS  PubMed  Google Scholar 

  35. Cocuzza S, Di Luca M, Maniaci A, Russo M, Di Mauro P, Migliore M, Serra A, Spinato G. Precision treatment of post pneumonectomy unilateral laryngeal paralysis due to cancer. FUTURE ONCOL. 2020;16(16s):45–53.

    CAS  PubMed  Google Scholar 

  36. Coussens LM, Zitvogel L, Palucka AK. Neutralizing tumor-promoting chronic inflammation: a magic bullet? Science. 2013;339(6117):286–91.

    CAS  PubMed  PubMed Central  Google Scholar 

  37. Varricchi G, Loffredo S, Marone G, Modestino L, Fallahi P, Ferrari SM, de Paulis A, Antonelli A, Galdiero MR. The Immune Landscape of thyroid Cancer in the context of Immune Checkpoint Inhibition. INT J MOL SCI 2019, 20(16).

  38. French JD, Kotnis GR, Said S, Raeburn CD, McIntyre RJ, Klopper JP, Haugen BR. Programmed death-1 + T cells and regulatory T cells are enriched in tumor-involved lymph nodes and associated with aggressive features in papillary thyroid cancer. J Clin Endocrinol Metab. 2012;97(6):E934–43.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. Ryu HS, Park YS, Park HJ, Chung YR, Yom CK, Ahn SH, Park YJ, Park SH, Park SY. Expression of indoleamine 2,3-dioxygenase and infiltration of FOXP3 + regulatory T cells are associated with aggressive features of papillary thyroid microcarcinoma. THYROID. 2014;24(8):1232–40.

    CAS  PubMed  Google Scholar 

  40. Liu Y, Yun X, Gao M, Yu Y, Li X. Analysis of regulatory T cells frequency in peripheral blood and tumor tissues in papillary thyroid carcinoma with and without Hashimoto’s thyroiditis. CLIN TRANSL ONCOL. 2015;17(4):274–80.

    PubMed  Google Scholar 

  41. Galdiero MR, Varricchi G, Loffredo S, Bellevicine C, Lansione T, Ferrara AL, Iannone R, di Somma S, Borriello F, Clery E, et al. Potential involvement of neutrophils in human thyroid cancer. PLoS ONE. 2018;13(6):e199740.

    Google Scholar 

  42. Qing W, Fang WY, Ye L, Shen LY, Zhang XF, Fei XC, Chen X, Wang WQ, Li XY, Xiao JC, et al. Density of tumor-associated macrophages correlates with lymph node metastasis in papillary thyroid carcinoma. THYROID. 2012;22(9):905–10.

    CAS  PubMed  PubMed Central  Google Scholar 

  43. Tanaka K, Kurebayashi J, Sohda M, Nomura T, Prabhakar U, Yan L, Sonoo H. The expression of monocyte chemotactic protein-1 in papillary thyroid carcinoma is correlated with lymph node metastasis and tumor recurrence. THYROID. 2009;19(1):21–5.

    CAS  PubMed  Google Scholar 

  44. Cunha LL, Morari EC, Nonogaki S, Marcello MA, Soares FA, Vassallo J, Ward LS. Interleukin 10 expression is related to aggressiveness and poor prognosis of patients with thyroid cancer. Cancer Immunol Immunother. 2017;66(2):141–8.

    CAS  PubMed  Google Scholar 

  45. Cunha LL, Morari EC, Nonogaki S, Soares FA, Vassallo J, Ward LS. Foxp3 expression is associated with aggressiveness in differentiated thyroid carcinomas. Clin (Sao Paulo). 2012;67(5):483–8.

    Google Scholar 

  46. Cunha LL, Marcello MA, Vassallo J, Ward LS. Differentiated thyroid carcinomas and their B7H1 shield. FUTURE ONCOL. 2013;9(10):1417–9.

    CAS  PubMed  Google Scholar 

  47. Liotti F, Kumar N, Prevete N, Marotta M, Sorriento D, Ieranò C, Ronchi A, Marino FZ, Moretti S, Colella R, et al. PD-1 blockade delays tumor growth by inhibiting an intrinsic SHP2/Ras/MAPK signalling in thyroid cancer cells. J Exp Clin Cancer Res. 2021;40(1):22.

    CAS  PubMed  PubMed Central  Google Scholar 

  48. Aghajani M, Graham S, McCafferty C, Shaheed CA, Roberts T, DeSouza P, et al. Clinicopathologic and prognostic significance of programmed cell death ligand 1 expression in patients with non-medullary thyroid cancer: a systematic review and meta-analysis. Thyroid. 2018, 28(3):349–361.

  49. Iyer PC, Dadu R, Gule-Monroe M, Busaidy NL, Ferrarotto R, Habra MA, Zafereo M, Williams MD, Gunn GB, Grosu H, et al. Salvage pembrolizumab added to kinase inhibitor therapy for the treatment of anaplastic thyroid carcinoma. J IMMUNOTHER CANCER. 2018;6(1):68.

    PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by the National Natural Science Foundation of China (No. 82271776).

Author information

Authors and Affiliations

  1. Department of Endocrinology, The First Affiliated Hospital, Sun Yat-sen University, Guangzhou, 510080, China

    Luyao Zhang, Shuang Yu, Shubin Hong, Xi Xiao, Zhihong Liao, Yanbing Li & Haipeng Xiao

Authors
  1. Luyao Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Shuang Yu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Shubin Hong
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xi Xiao
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Zhihong Liao
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yanbing Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Haipeng Xiao
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

All authors contributed to the study conception and design. HX designed the research flow. LZ and SY did the data collection and validation. LZ, XX and SH analyzed the data. LZ drafted the manuscript. ZL and YL revised the manuscript critically. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Haipeng Xiao.

Ethics declarations

Ethics approval

This study was performed in line with the principles of the Declaration of Helsinki. Approval was granted by the First Affiliated Hospital of Sun Yat-sen University ([2021]109).

Consent for publication

Not applicable.

Conflict of interest

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Supplementary Material 3

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, L., Yu, S., Hong, S. et al. Comprehensive analysis of BTNL9 as a prognostic biomarker correlated with immune infiltrations in thyroid cancer. BMC Med Genomics 16, 234 (2023). https://doi.org/10.1186/s12920-023-01676-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-023-01676-8

Keywords

BMC Medical Genomics

ISSN: 1755-8794

Contact us