Prognostic biomarker TUBA1C is correlated to immune cell infiltration in the tumor microenvironment of lung adenocarcinoma

TUBA1C is a microtubule component that is involved in a variety of cancers. Our main objective was to investigate TUBA1C expression, its prognostic value, its potential biological functions, and its impact on the immune system of patients with lung adenocarcinoma (LUAD). The Cancer Genome Atlas (TCGA), Gene Expression Profiling Interactive Analysis (GEPIA) and Immunohistochemistry Analysis were used to analyze TUBA1C expression, its clinicopathology, overall survival (OS), and disease-free survival (DFS) in LUAD patients. We also determined the correlation between TUBA1C and tumor-infiltrating immune cells (TIICs) by using CIBERSORT and GEPIA databases. To determine the expression of TUBA1C in LUAD, we analyzed a collection of immune infiltration levels and cumulative survival of LUAD tissues in TIMER database. By using UALCAN, STRING, and GeneMANIA databases, we investigated the protein-coding genes related to TUBA1C and its co-expression genes in LUAD tissues. Gene set enrichment analysis (GSEA) was performed by using the TCGA dataset. The mRNA and the protein expression of TUBA1C were found to be up-regulated in LUAD tissues. The univariate analysis indicated that an increased expression of TUBA1C was significantly correlated to the following parameters: age, stage, and lymph node metastasis. An over-expression of TUBA1C was associated with a poor prognosis of LUAD. In TIMER and CIBERSORT databases, we found that TUBA1C is correlated with 13 types of TIICs: activated B cell, activated CD4 T cell, central memory CD4 T cell, effector memory CD8 T cell, eosinophils, immature B cell, gamma-delta T cell, immature dendritic cell, mast cell, memory B cell, natural killer T cell, regulatory T cell, and type 2T helper cell. By performing GSEA, we found that TUBA1C is closely correlated to cell cycle, p53 signaling pathway, glycolysis, and gluconeogenesis. Our findings indicate that TUBA1C is associated with TIICs in tumor microenvironment. Therefore, it serves as a novel prognostic biomarker and a target for future treatment methods of LUAD.


Background
The morbidity and mortality of lung cancer has been the highest among all types of cancers, so it is the most common form of malignancy worldwide [1]. Lung adenocarcinoma (LUAD) is a major pathological subtype of non-small cell lung cancer (NSCLC), which approximately accounts for 40% of lung cancer cases [2]. Surgery, chemotherapy, radiotherapy, and targeted therapy are the conventional methods of treating patients with NSCLC; however, the prognosis of these patients has been unsatisfactory till date [3]. Since the past decade, immune checkpoint inhibitors (ICIS) have been used to treat patients with NSCLC, and they have modified the treatment pattern of this refractory disease [4]. Tumor-infiltrating immune cells (TIICs) have impacted the immune system and tackled abnormal biological behaviors in a complex way, so they play a key role in eliciting the body's response to immunotherapy [5].
Microtubule is an important component of the eukaryotic cytoskeleton. Moreover, it is one of the most functional proteins, which play an important role in dynamic polymerization and depolymerization through cell replication and division [6]. Microtubules are uniformly assembled from highly conserved α/β-tubulin heterodimers [7]. Recently, several studies have reported that α-tubulin is involved in the occurrence of a variety of tumors, such as lung cancer, breast cancer and prostatic cancer [8][9][10]. In addition, α-tubulin is also associated with the development of astrocytoma and chemotherapeutic resistance of liver cancer [11]. Moreover, TUBA1C is a subtype of α-tubulin, and its overexpression is associated with the poor prognosis of hepatocellular carcinoma (HCC) and pancreatic ductal adenocarcinoma [12,13]. However, no previous study has elucidated how the overexpression of TUBA1C affects LUAD patients.
In this study, IHC and Gene Expression Profiling Interactive Analysis (GEPIA) were performed to elucidate the correlation between TUBA1C and LUAD. Furthermore, we used the computational algorithms CIBERSORT and TIMER to explore the relationship between TUBA1C and TIICs in LUAD patients. In addition, STRING software, GeneMANIA Analysis, and Gene Set Enrichment Analysis (GSEA) were used to further study the function and mechanism of TUBA1C in LUAD patients. The findings proved that TUBA1C played an important role in the development of LUAD. Furthermore, we elucidated the potential relationship between TUBA1C and tumorimmune interaction.

Data acquisition
The TCGA database (https ://porta l.gdc.cance r.gov/) was used to obtain the data of TUBA1C expression in LUAD tissues and normal tissues. A total of 535 tumor tissues (Additional file 1) and 59 normal tissues were included in the data analysis. Meanwhile, 390 LUAD tissues and 60 normal tissues were treated from the Affiliated Hospital of Nantong University in China. All patients were treated by surgical resection between 2012 and 2013. All clinical data on the patients were carefully recorded after the diagnosis of LUAD by two pathologists. The follow up was 60 months. All experiments involving patient specimens were approved by the Ethics Committee of the Affiliated Hospital of Nantong University, China.
To determine the expression of TUBA1C in LUAD tissues, we performed GEPIA and constructed boxplot, pathological stage plot, and survival curves, such as overall survival (OS) curve and disease-free survival (DFS) curve. Meanwhile, the relationship between TUBA1C and immune infiltrating cells was analyzed by GEPIA. The Spearman method was used to determine the correlation coefficient of the relationship.

Immunohistochemistry analysis
An immunohistochemistry (IHC) assay was conducted as previously described [15]. Briefly, the LUAD samples were deparaffinized and rehydrated. The primary antibody was that against TUBA1C (1:100 dilution; ab222849; abcam). The positive expression of TUBA1C was localized in cytoplasm. The scoring criteria for IHC staining were based on the intensity of the stain and the percentage of immunoreactive cells, as previously described [15].

TIMER analysis
In this experiment, we used TIMER (https ://cistr ome. shiny apps.io/timer /) algorithm to comprehensively elucidate the correlation between different types of tumors and TIICs [16]. Moreover, we performed a series of analyses to determine the expression of TUBA1C in different types of cancer and to comprehend its correlation with the abundance of TIICs. These TIICs were further classified as follows: B cells, CD4+ T cells, CD8+ T cells, macrophages, neutrophils, and dendritic cells. By using the TIMER algorithm, we performed a correlation analysis between TUBA1C, its related genes, and the markers of immune cells [17]. The expression of TUBA1C gene was plotted on the x-axis, and the expression of its related marker genes was plotted on the y-axis.

CIBERSORT analysis
The computational method of CIBERSORT [18] (http:// ciber sort.stanf ord.edu/) is a deconvolution algorithm based on gene expression, and it is used to evaluate the changes of a group of genes relative to all other genes in a sample. Based on the expression of TUBA1C, we classified 497 samples into high and low expression groups. Using CIBERSORT algorithm, we measured the immune responses of 28 TIICs and evaluated the relationship between these TIICs and the expression of TUBA1C in LUAD tissues. Our main goal was to determine the correlation between these TIICs.

UALCAN analysis
The University of Alabama Cancer Database (UALCAN) (http://ualca n.path.uab.edu/index .html) is a visual portal of the Cancer Genome Atlas (TCGA) database [19]. This database was used to analyze the positively and negatively expressed protein-coding genes related to TUBA1C in LUAD tissues.

STRING and GeneMANIA analysis
To obtain the information on protein-protein interaction, most scientists prefer using the following two datasets: STRING and GeneMANIA [20,21]. The protein-protein interaction network of TUBA1C was predicted with the help of STRING and GeneMANIA datasets.

Gene set enrichment analysis (GSEA)
Gene set enrichment analysis (GSEA) was performed on the normalized RNA-Seq data, which was obtained from TCGA database [22]. The gene ontology (GO) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were used to investigate the possible biological functions of TUBA1C. A false discovery rate (FDR) < 0.050 and a nominal P < 0.050 were considered to be statistically significant.

Statistical analysis
All statistical analyses were performed by using R language (version 3.5.3). Kaplan-Meier univariate analysis was used to determine the effect of TUBA1C expression on survival. Multivariate Cox analysis was performed to determine the expression of TUBA1C and the effect of other pathological and clinical factors (age, gender, stage, tumor status, and lymph node) on overall survival (OS). The results were considered to be statistically significant when P < 0.05.

The mRNA expression of TUBA1C in different tumors
In order to determine the difference between the expression of TUBA1C in tumor tissues and normal tissues, we used TIMER algorithm to analyze the mRNA expression of TUBA1C in different types of tumors that were obtained from TCGA database (Fig. 1a). Compared to normal tissues, the mRNA expression of TUBA1C was found to be higher in the following types of tumors: bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma

The levels of TUBA1C predicts the prognosis of LUAD
The TCGA database was used to identify the differences in the mRNA levels of TUB1AC in normal tissues and LUAD tissues. The data obtained from 535 tumor tissues and 59 normal tissues was analyzed in this study. As shown in Fig. 1b, a boxplot was plotted using the data of expression of TUBA1C in normal and tumor tissues.
The results indicate that the expression of TUBA1C in tumor tissues was significantly higher than in normal tissues (P = 8.913e-21). As shown in Fig. 1c, the Kaplan-Meier (KM) survival curve showed a high expression of TUBA1C in LUAD tissues, which indicated poor prognosis of LUAD. As shown in Fig. 1d, multivariable analysis was performed by adjusting for the following factors: age, stage, lymph node metastasis, and TUBA1C expression. Furthermore, boxplot, pathological stage plot, and survival curves were constructed by plotting the expression levels of TUBA1C with the help of GEPIA. As shown in Fig. 2, a high expression of TUBA1C was significantly associated with the following parameters: different disease states (Tumor or Normal) (P < 0.05), pathological stage (P = 1.82e-04), overall survival (P = 3.8e−05), and disease-free survival (P = 0.047). Immunohistochemical analysis of TUBA1C expression in 390 cases of lung adenocarcinoma and 60 cases of normal tissues is shown in  Table 1). The KM curve analysis of TUBA1C shows the same results as TCGA data (Fig. 3g). Multivariable analysis demonstrated that only the TUBA1C expression was an independent prognostic factor in lung adenocarcinoma (Fig. 3h).

The relationship between the expression of TUBA1C and TIICs
Several evidences prove that the characteristics of TIICs were significantly correlated to the occurrence and development of tumor tissues [23][24][25]. Therefore, we explored whether the expression of TUBA1C was associated with TIICs in LUAD tissues. Based on the expression of TUBA1C, we segregated 497 samples into high and low expression groups. A total of 248 high expression samples and 249 low expression samples met the screening criteria. The two groups exhibited differences in the 28 proportions of immune cells, which were observed by downloading the gene expression source from an established and trusted computing resource (CIBERSORT, Newman et al., Stanford University, USA). Figure 4a displays the results of the 28 immune cell subsets. The expression of TUB1AC was significantly correlated to the following types of immune cells: activated B cell, activated CD4 T cell, central memory CD4 T cell, effector memory CD8 T cell, eosinophils, gamma-delta T cell, immature B cell, immature dendritic cell, mast cell, memory B cell, natural killer T cell, regulatory T cell, and type 2 T-helper cell. In particular, activated CD4 T cell, effector memory CD8 T cell, gamma-delta T cell, memory B cell, natural killer T cell, and regulatory T cell were present in higher proportions in the high expression group than in others. Moreover, activated B cell, central memory CD4 T cell, eosinophils, immature B cell, mast cell, and type 2 T -helper cell were also present in higher proportion in the high expression group (Fig. 4a). Meanwhile, we used TIMER and GEPIA databases to further investigate the relationship between TUBA1C and the diverse set of immune infiltrating cells. After adjusting the correlation coefficients according to purity, we found that the expression of TUBA1C was significantly correlated to T cells, B cells, natural killer cells, and neutrophils in LUAD tissues (Table 2). Furthermore, we used the GEPIA database to analyze the correlation between TUBA1C expression and the above-mentioned  (Table 3).

In the TIMER database, TUBA1C expression was associated with immune infiltration level in LUAD and cumulative survival in LUAD
In this study, TIICs played a decisive role in the prognosis and survival of LUAD patients [26]. Therefore, we used TIMER database to further explore the relationship between the prognosis and survival of infiltrating immune cells and TUBA1C expression in LUAD tissues. As shown in Fig. 4a, b negative correlation exists between TUBA1C expression levels and the infiltrating level of B cells (r = − 0.331, P = 8.22e−14) and CD4+ T cells (r = − 0.2, P = 9.64e−06). In LUAD tissues, TUBA1C expression was associated with poor prognosis and high immune infiltration. As shown in Fig. 4c, the infiltrating level of B cell and TUBA1C expression were the factors related to the cumulative survival rate of LUAD.

The protein-coding genes of TUBA1C and its co-expression genes in LUAD
To investigate the potential molecular mechanisms through which TUBA1C elicits tumorigenesis in LUAD, we identified the protein-coding genes of TUBA1C and its co-expression genes in LUAD. As shown in Fig. 5, UALCAN database was used to identify genes that showed a positive and negative correlation with TUBA1C in LUAD tissues. The two-sided Pearson's correlation coefficient analysis and z-test were performed by using R language, which is based on the gene expression data extracted from TCGA. The top 10 protein-coding genes that positively correlated with TUBA1C were as follows: TUBA1B, PLK1, BIRC5, CCNA2, CCNB1, EPR1, MAD2L1, RAN, SKA3, and PACGAP1. On the other hand, the top ten protein-coding genes that negatively correlated with TUBA1C were as follows: CIRBP, VAMP2, CBX7, NICN1, C5orf53, FRY, CRY2, GNG7, CD302, and FXYD1. Furthermore, STRING and Gene-MANIA tools were used to analyze the interaction between TUBA1C and protein-coding genes mentioned earlier. Figure 6 illustrates the results of the analysis.

Gene sets enriched in TUBA1C expression
In this experiment, GO and KEGG pathway analysis were performed to determine the potential biological functions of TUBA1C. We observed a significant difference in the enrichment of GO and KEGG pathways, which were used to analyze samples that showed a high expression of TUBA1C (FDR < 0.050, P < 0.050). Table 4 presents 10 KEGG and GO pathways that were associated with a high expression of TUBA1C. As shown in Fig. 7a, the 10 KEGG pathways that were positively correlated to the high expression of TUBA1C were as follows: cell cycle, p53 signaling pathway, basal transcription factors, ubiquitin mediated proteolysis, glycolysis gluconeogenesis, citrate cycle TCA cycle, oxidative phosphorylation, pancreatic cancer, renal cell carcinoma, bladder cancer. As shown in Fig. 7b, GO analysis also revealed the following ten positively correlated categories: cell cycle G2/M phase transition, cellular response to oxygen levels, chromosome segregation, interleukin 1 mediated signaling pathway, multi-organism localization, condensed chromosome, ubiquitin like protein conjugating enzyme activity, microtubule cytoskeleton organization involved in mitosis, sister chromatid segregation, and mitotic spindle organization. These results indicate that the pathways of cell cycle, glycolysis, and gluconeogenesis are involved in the carcinogenesis of LUAD, which is strongly related to TUBA1C expression.

Discussion
In molecular biology, TUBA1C is a kind of α-tubulin subtype that is related to microtubules. It is a multifunctional cytoskeleton protein, and it participates in the process of cell mitosis and cell division [27,28]. Previous studies have reported that when the expression of TUBA1C is up-regulated, it significantly affects the growth and progression of tumor cells. This indicates that TUBA1C plays a pivotal role in the proliferation and cell cycle of various tumors [29,30]. Recent studies have investigated the structure and function of microtubules (MT). These studies have reported its potential role in innate and adaptive immune systems [31]. Moreover, TIICs are found to be involved in the growth, invasion, and metastasis of lung cancer [32]. However, no previous study has reported about the relationship between TUBA1C expression and TIICs in LUAD tissues.
In this study, we investigated TUBA1C expression levels in various types of tumors, which were identified from the TCGA database by using TIMER algorithm.
Compared to normal tissues, TUBA1C expression was found to be higher in the following types of tumors: BLCA, BRCA, CHOL, COAD, ESCA, HWSC, KIRC, KIRP, LIHC, LUAD, LUSC, PRAD, READ, STAD, THCA, and UCEC (Fig. 1a). The data of LUAD patients was acquired from TCGA database and used to estimate the prognosticative value of TUBA1C. Our main goal was to figure out whether TUBA1C can be used as a prognosticative biomarker. Furthermore, we found that the expression of TUBA1C was extremely higher in LUAD than in normal tissues. (Fig. 1b). As shown in Figs. 1c and 3g, the KM survival curve exhibited a high expression of TUBA1C in LUAD tissues, which indicated a poor prognosis of LUAD. Age, tumor grade, lymph node metastasis, and TUBA1C expression were considered as independent prognostic factors in multivariable analysis, which is illustrated as a forest boxplot in Fig. 1d. By using the GEPIA database, we found that TUBA1C expression was significantly related to the following parameters: different disease states (tumor or normal) (P < 0.05), pathological stage (P = 1.82e−04), overall survival (P = 3.8e−0.5), and disease-free survival (P = 0.047) (Fig. 2). By performing immunohistochemical analysis, we found that the Fig.4 The relationship between the expression of TUBA1C and TIICs. a The results of the relative ratios of TIIC that were obtained using the CIBERSORT algorithm. The ratio of 28 immune cells in LUAD tissues in the TUBA1C high and low expression groups (*P < 0.05, **P < 0.01, ***P < 0.001). b Negative correlation exists between the TUBA1C expression level and infiltrating levels of B cell (r = − 0.331, P = 8.22e−14) and CD4+ T cells (r = − 0.2, P = 9.64e−06) in LUAD. c Cumulative survival is related to B cell and the expression of TUBA1C in LUAD. (The B cell and the expression of TUBA1C are factors related to the cumulative survival rate of LUAD over time) expression level of TUBA1C was higher in tumor tissues than in normal tissues (Fig. 3a-f ). From our data, there is no correlation between the level of TUBA1C expression and clinicopathological parameters of lung adenocarcinoma (Table 1). Multivariable analysis demonstrated that only the TUBA1C expression was an independent prognostic factor in lung adenocarcinoma (Fig. 3h). All these findings suggest that TUBA1C is a prognostic biomarker of LUAD.
In this study, we found that various immune infiltration levels of LUAD could induce different expression levels of TUBA1C. By using CIBERSORT algorithm, we found that TUBA1C expression is strongly associated with infiltration levels of following immune cells: activated B cell, activated CD4 T cell, central memory CD4 T cell, effector memory CD8 T cell, eosinophils, gamma-delta T cell, immature B cell, immature dendritic cell, mast cell, memory B cell, natural killer T cell, regulatory T cell and type 2 T-helper cell (Fig. 4a). Furthermore, TIMER and GEPIA databases showed how gene markers of different immune cells were related to TUBA1C expression. This may indicate that the tumor immune microenvironment (Tables 2, 3) of LUAD was regulated by TUBAC1. It is a well-known fact that T cells and B cells comprise about two-thirds of lung TIICs, while the remaining lung TIICs are composed of tumor-associated macrophages (TAMs) and a small number of infiltrating dendritic cells and natural killer cells (NK) [33]. Previous studies have thoroughly investigated the role of tumor-associated T cells in the development of lung cancer. These studies have reported that CD4 + Th1 cells and activated CD8 + T cells often elicited type I immune responses, which indicate a favorable prognosis of LUAD [34,35]; however, Th2, Th17, and Foxp3 + regulatory T (Treg) cells were found to be associated with tumor progression   [37][38][39].
Many studies have proved that the expression level of many genes is related to poor prognosis and TIICs. By analyzing the TIMER database, we found that the expression of TUBA1C is negatively correlated to B  (Fig. 4b, c). Although TUBA1C is a tubulin, we still do not know the mechanism through it regulates LUAD. First, we analyzed the protein-coding genes related to TUBA1C and its co-expression genes in LADC tissues. The top 10 protein-coding genes that positively correlated with TUBA1C are as follows: TUBA1B, PLK1, BIRC5, CCNA2, CCNB1, EPR1, MAD2L1, RAN, SKA3, and PACGAP1. On the other hand, the top 10 negatively correlated genes are as follows: CIRBP, VAMP2, CBX7, NICN1, C5orf53, FRY, CRY2, GNG7, CD302, and FXYD1 (Fig. 5). Furthermore, STRING and Gene MANIA databases illustrated the protein interaction between TUBA1C and other partners (Fig. 6). The proteins related to TUBA1C perform following biological functions: they regulate the cell cycle, mitosis, DNA damage response, cell proliferation, and aging. Thereafter, GO and KEGG pathway analysis revealed that an up-regulated expression of TUBA1C is primarily related to cell cycle, p53 signaling pathway, glycolysis, and gluconeogenesis (Fig. 7, Table 4). Previous studies have also reported that TUBA1C is associated with cell proliferation, and it also regulates cell cycle in many types of cancers. Furthermore, the expression of TUBA1C was reported to be correlated with p53 expression in pancreatic ductal adenocarcinoma [12]. Previous studies have reported that tubulin regulates cell metabolism and glucose stress response, reducing the dependence of cells on glycolysis. These events ultimately promote cell survival in cancer patients [40]. These results are helpful to understand the biological role played by TUBA1C in the development of LUAD. In future clinical practice, the expression of TUBA1C in lung adenocarcinoma tissue may be used to predict the prognosis and the efficacy of immunotherapy of patients.

Conclusion
In conclusion, this is the first report to prove that TUBA1C is a new marker of LUAD. This work further proved that TUBA1C played a pivotal role in the cell cycle and immune microenvironment of LUAD. By further understanding its range of functions, we can make