DNAJB11 predicts a poor prognosis and is associated with immune infiltration in thyroid carcinoma: a bioinformatics analysis

Abstract

Objective

To investigate the prognostic value of the co-chaperone protein DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11) in thyroid carcinoma (THCA).

Methods

This bioinformatics analysis study evaluated the prognostic value of DNAJB11 mRNA levels in THCA based on data from The Cancer Genome Atlas (TCGA). The levels of DNAJB11 mRNA in THCA and normal tissues were compared with Wilcoxon signed rank test. Kaplan–Meier survival curve analysis and Cox regression analysis were performed to evaluate the correlation between DNAJB11 mRNA levels and survival. Gene Ontology (GO) enrichment analysis was used to elucidate the functional enrichment difference.

Results

Data from the 502 patients with THCA from the TCGA database were analysed. DNAJB11 mRNA was downregulated in THCA tissues compared with normal tissues. Decreased levels of DNAJB11 mRNA were significantly correlated with T stage, N stage, pathological stage, histological type, extrathyroidal extension and BRAF gene status. The low levels of DNAJB11 mRNA were associated with a shorter progression-free interval. GO enrichment analysis showed that DNAJB11 was involved in immune-related biological processes.

Conclusion

Low levels of DNAJB11 mRNA were associated with poor prognosis in THCA.

Keywords

DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11)thyroid carcinoma prognosis The Cancer Genome Atlas Gene Set Enrichment Analysis

Introduction

Thyroid carcinoma (THCA) is the most common endocrine malignant tumour, with an increasing global incidence each year.¹ THCA is classified into four subtypes based on histology: papillary thyroid carcinoma, follicular thyroid carcinoma, medullary thyroid carcinoma and anaplastic thyroid carcinoma.² Papillary thyroid carcinoma is responsible for 80% of thyroid carcinomas, which are more common in young women.³ Despite advances in detection and surgical management, such as surgical resection, radiotherapy and levothyroxine treatment, metastasis and recurrence are common, with a 5-year survival rate of approximately 59% in patients at an advanced stage.⁴ As a result, it is critical to identify reliable predictors of prognostic stratification and provide new therapeutic targets. Many different types of molecular signals (such as DNA mutations and abnormal mRNA and microRNA expression) have been identified independently as favourable or unfavourable prognostic biomarkers.⁵

The DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11) gene is a protein coding gene that encodes a co-chaperone protein involved in protein processing in the endoplasmic reticulum (ER).⁶ The protein performs multifaceted functions involved in coordinating ER and extracellular proteostasis.^7,8 DNAJB11 binds to misfolded proteins in the ER and transports them to the ER HSP70 Bip for ATP-dependent chaperoning.^9–11 Recent research has revealed that the endoplasmic reticulum co-chaperone ERdj3/DNAJB11 promotes hepatocellular carcinoma progression through suppressing AATZ degradation.¹² Moreover, DNAJB11 is positively correlated with cirrhosis, advanced hepatocellular carcinoma status and poor survival outcomes.¹² DNAJB11 is found in high concentrations in the cytoplasm of oral cavity squamous cell carcinoma cells.¹³ The ER stress response has been implicated in malignant cell survival and development.¹⁴ DNAJB11 protects against ER stress by stabilizing proteins.¹⁵ DNAJB11 is associated with a variety of signalling pathways, including aberrant signalling pathways associated with cancer.¹⁶

This current study used bioinformatics analysis of publicly available sequencing data to gain a better understanding of the potential actions of DNAJB11 in THCA and its regulatory network. RNA sequence data from The Cancer Genome Atlas (TCGA) were used in this current study to assess the prognostic value of DNAJB11 gene expression in human THCA. The differences in the levels DNAJB11 mRNA between tumours from patients with THCA and normal tissues were analysed so that any correlation between DNAJB11 mRNA levels and clinicopathological parameters could be determined.

Materials and methods

Data acquisition and preprocessing

The RNA sequence data and the corresponding clinicopathological information of patients with THCA were obtained from the TCGA database (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) to forecast DNAJB11 mRNA levels; and the association between the levels and the overall survival (OS) of patients with THCA. The number of fragments per kilobase of non-overlapped exons per million fragments mapped (FPKM) was calculated and then the FPKM was subsequently transferred into transcripts per kilobase million (TPM) values. There were 502 patients with THCAs included in the further analysis, while patients without RNA sequence data (n = 10) and those with OS <30 days (n = 20) were excluded. As this study analysed bioinformatics data from a public database and did not involve human subjects, informed consent was not required. In addition, no approval from an institutional review board was required.

Differentially expressed gene (DEG) analysis

The DEGs between THCA and normal tissues were identified using HTSeq-Counts by the DESeq2 software by using R.¹⁷ The Benjamini–Hochberg procedure was introduced to reduce the false positive rate (FDR) in multiple comparisons. The cut-off criteria were set as |log2 Fold Change (FC)| > 1.5 and adjusted P-value <0.05.

GO enrichment analysis

Metascape online software (http://metascape.org)¹⁸ was used as a tool to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DNAJB11-related DEGs. The significantly enriched (P-value <0.01, minimum count of 3, and enrichment factor of >1.5) terms of biological process were visualized.

GSEA enrichment analysis

Gene Set Enrichment Analysis (GSEA) was performed using the gseGO and gseKEGG functions of R package clusterProfiler (3.8.0).¹⁹ In this study, GSEA was used to elucidate the significant function and pathway differences between high- and low-DNAJB11 groups. The analysis process was repeated 1000 times using the default weighted enrichment statistics methods. Genes with an adjusted P-value <0.05 and FDR q-value <0.25 were considered significantly related genes.

Immune infiltration analysis by ssGSEA

The immune infiltration analysis of DNAJB11 was done by applying a single-sample Gene Set Enrichment Analysis (ssGSEA) method from the Gene Set Variation Analysis (GSVA) package²⁰ in R (version 3.6.3) for 24 types of immune cells in tumour samples. The following 24 types of immune cells were obtained: CD4+ T cells, CD8+ T cells, gamma delta T cells, B cells, CD56bright natural killer cell, CD56dim natural killer cells, natural killer cells, dendritic cells (DCs), activated dendritic cells (aDCs), immature dendritic cells (iDCs), plasmacytoid dendritic cells, T helper cells, regulatory T (T_reg) cells, follicular helper T cells, type 1 helper T (T_h1) cells, type 2 helper T (T_h2) cells, type 17 helper T cells, central memory T (Tcm) cells, effector memory T cells (Tem), macrophages, mast cells, neutrophils, eosinophils and cytotoxic cells.

Statistical analyses

All statistical analysis and plots were analysed using the R statistical package (R version 3.6.3; R Foundation for Statistical Computing, Vienna, Austria). The relationship between clinicopathological features and the levels of DNAJB11 were analysed using Wilcoxon signed rank test. The Kaplan–Meier method and Cox regression analysis were used to evaluate prognostic factors. Multivariate Cox analysis was used to compare the influence of DNAJB11 levels on survival along with other clinical factors, including stage, grade, histological type, extrathyroidal extension, v-raf murine sarcoma viral oncogene homolog B1 (BRAF) gene status, neuroblastoma RAS viral oncogene homolog (NRAS) gene status and Harvey rat sarcoma viral oncogene homolog (HRAS) gene status. The median DNAJB11 level was regarded as the cut-off value. The hazard ratio (HR) with 95% confidence interval (CI) was measured to estimate the hazard risk of individual prognostic factors. A P-value <0.05 was considered statistically significant in all tests.

A nomogram was constructed based on the results of the multivariate analysis to predict the prognosis of patients with THCA by R packages rms. Multivariate Cox regression analysis and Akaike’s information criterion method were used to evaluate each variable. The TCGA THCA cohort was stratified into a high-risk group and low-risk group based on the median risk score. The receiver operating characteristic (ROC) curve was used to evaluate the predictive accuracy of the prognostic model.

Results

The data collected from the TCGA included 502 patients with THCA that had both clinical and gene expression data. The characteristics of the patients are shown in Table 1. The cancer status included 502 tumour tissues and 58 adjacent tissues. A total of 367 female patients and 135 male patients were analysed in this study. The median age at diagnosis was 46 years. The number of patients with pathological stage I, II, III and IV were 280 (55.8%), 55 (11.0%), 112 (22.4%) and 55 (11.0%). A total of 223 of 502 (44.4%) patients had lymph node metastases and nine of 502 (1.8%) had distant metastases.

Table 1.

Clinical and demographic characteristics of the patients (n = 502) with thyroid carcinoma for whom data were collected from The Cancer Genome Atlas (TCGA) database stratified according to the levels of DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11) mRNA.

Characteristic	Level	n ^#	Low levels of DNAJB11 mRNA	High levels of DNAJB11 mRNA	Statistical analysis^a	Test
T stage	T1	143	62 (43.4%)	81 (56.6%)	P < 0.001	Exact
	T2	164	63 (38.4%)	101 (61.6%)
	T3	170	105 (61.8%)	65 (38.2%)
	T4	25	19 (76.0%)	6 (24.0%)
N stage	N0	229	97 (42.4%)	132 (57.6%)	P < 0.001
N stage	N1	223	132 (59.2%)	91 (40.8%)
M stage	M0	282	153 (54.3%)	129 (45.7%)	P = 0.01	Exact
M stage	M1	9	5 (55.6%)	4 (44.4%)
Pathological stage	Stage I	280	120 (42.9%)	160 (57.1%)	P < 0.001
	Stage II	55	21 (38.2%)	34 (61.8%)
	Stage III	112	73 (65.2%)	39 (34.8%)
	Stage IV	55	35 (63.6%)	20 (36.4%)
Residual tumour	R0	384	185 (48.2%)	199 (51.8%)	P = 0.03	Exact
	R1	52	34 (65.4%)	18 (34.6%)
	R2	4	3 (75.0%)	1 (25.0%)
Histological type	Classical	356	182 (51.1%)	174 (48.9%)	P < 0.001	Exact
	Follicular	101	29 (28.7%)	72 (71.3%)
	Other	9	6 (66.7%)	3 (33.3%)
	Tall cell	36	34 (94.4%)	2 (5.6%)
Sex	Female	367	191 (52.0%)	176 (48.0%)	NS
Sex	Male	135	60 (44.4%)	75 (55.6%)
Race	Asian	51	21 (41.2%)	30 (58.8%)	NS
	Black or African American	27	17 (63.0%)	10 (37.0%)
	White	332	179 (53.9%)	153 (46.1%)
Neoplasm location	Bilateral	86	42 (48.8%)	44 (51.2%)	NS
	Isthmus	22	9 (40.9%)	13 (59.1%)
	Left lobe	175	82 (46.9%)	93 (53.1%)
	Right lobe	213	114 (53.5%)	99 (46.5%)
Primary neoplasm focus type	Multifocal	226	107 (47.3%)	119 (52.7%)	NS
Primary neoplasm focus type	Unifocal	266	141 (53.0%)	125 (47.0%)
Thyroid gland disorder history	Lymphocytic thyroiditis	71	33 (46.5%)	38 (53.5%)	NS
	Nodular hyperplasia	68	30 (44.1%)	38 (55.9%)
	Normal	280	143 (51.1%)	137 (48.9%)
	Other, specify	25	15 (60.0%)	10 (40.0%)
Extrathyroidal extension	No	331	138 (41.7%)	193 (58.3%)	P < 0.001
Extrathyroidal extension	Yes	153	105 (68.6%)	48 (31.4%)
BRAF status	Mut	285	191 (67.0%)	94 (33.0%)	P < 0.001
BRAF status	WT	200	53 (26.5%)	147 (73.5%)
NRAS status	Mut	39	8 (20.5%)	31 (79.5%)	P < 0.001
NRAS status	WT	446	236 (52.9%)	210 (47.1%)
HRAS status	Mut	17	1 (5.9%)	16 (94.1%)	P < 0.001	Exact
HRAS status	WT	468	243 (51.9%)	225 (48.1%)
Age, years			50.00 (38.00–60.00)	43.00 (33.00–56.00)	P = 0.001	Non-norm

Data presented as n of patients (%) or median (interquartile range).

^aWilcoxon signed rank test.

^#Data are missing for some patients.

T, tumour; N, node; M, metastasis; BRAF, v-raf murine sarcoma viral oncogene homolog B1; MUT, mutant; WT, wild type; NRAS, neuroblastoma RAS viral oncogene homolog; HRAS, Harvey rat sarcoma viral oncogene homolog; NS, no significant between-group difference (P ≥ 0.05).

There were significantly lower levels of DNAJB11 mRNA in THCA tissues compared with non-paired normal tissues and paired adjacent samples (P < 0.001 for both comparisons) (Figures 1A and 1B). The ROC curve analysis showed that DNAJB11 mRNA levels can effectively discriminate THCA from normal tissue in the TCGA database. The area under the curve of DNAJB11 mRNA level was 0.739 (Figure 1C). Compared with the normal tissues, differential expression was observed in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), oesophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), sarcoma (SARC), stomach adenocarcinoma (STAD), THCA and uterine corpus endometrial carcinoma (UCEC) (P < 0.05 for all comparisons) (Figure 1D).

Figure 1.

Identification of the differentially expressed DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11) gene. (a) Levels of DNAJB11 mRNA from The Cancer Genome Atlas (TCGA) database in non-paired thyroid carcinoma (THCA) and normal samples were analysed using Wilcoxon signed rank test. Lower levels of DNAJB11 mRNA were observed in THCA compared with normal samples. (b) Levels of DNAJB11 mRNA from the TCGA database in paired THCA and adjacent samples were analysed using Wilcoxon signed rank test. Lower levels of DNAJB11 mRNA were observed in THCA compared with adjacent samples. (c) Receiver operating characteristic curve analysis of DNAJB11 mRNA levels showing promising discriminatory power between THCA and non-tumour tissues. (D) Human DNAJB11 levels in different tumour types from the TCGA database. Wilcoxon signed rank test, ns, P ≥ 0.05; *P < 0.05; **P < 0.01; ***P < 0.001. ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; DLBC, diffuse Large B-cell; ESCA, oesophageal 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 leukaemia; LGG, low grade glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma; OV, ovarian cancer; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PRAD, prostate adenocarcinoma; READ, rectum adenocarcinoma; SARC, sarcoma; SKCM, ▪; STAD, stomach adenocarcinoma; TGCT, testicular germ cell tumour; THYM, thymoma; UCEC, uterine corpus endometrial carcinoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma. The colour version of this figure is available at: http://imr.sagepub.com.

The Kruskal-Wallis test and Wilcoxon signed rank test showed that lower levels of DNAJB11 mRNA were significantly correlated with T stage (P < 0.001), N stage (P < 0.001), pathological stage (P < 0.001), histological type (P < 0.001), extrathyroidal extension (P < 0.001) and BRAF status (P < 0.001) (Figure 2).

Figure 2.

Correlation between the levels of DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11) mRNA and clinicopathological characteristics: (a) tumour (T) stage; (b) node (N) stage; (c) pathological stage; (d) histological type; (e) extrathyroidal extension; (f) v-raf murine sarcoma viral oncogene homolog B1 (BRAF) gene status. Wilcoxon signed rank test. The heavy central black horizontal lines for each group are the means. The extremities of the box are the 25th and 75th percentiles. The error bars represent the minimum and maximum outliers. The colour version of this figure is available at: http://imr.sagepub.com.

To explore the prognostic value of DNAJB11 mRNA levels in THCA outcome, the median value of DNAJB11 mRNA was used to divide patients into DNAJB11-high and DNAJB11-low level groups. The survminer R package was used to perform a Kaplan–Meier analysis.²¹ Compared with the DNAJB11-high level group, the DNAJB11-low level group was significantly associated with shorter progression-free survival (PFS) (HR 0.41; 95% CI 0.22, 0.74; P = 0.003) (Figure 3A). To provide clinicians with a quantitative method for predicting the prognosis of THCA patients, a nomogram was constructed to integrate DNAJB11 mRNA level and other independent prognostic variables including pathological stage and age (Figure 3B). The C-index for the nomogram was 0.678 with 1000 bootstrap replicates (95% CI 0.644, 0.713). The bias-corrected line in the calibration plot was found to be close to the ideal curve (the 45-degree line), which indicated good agreement between the prediction and the observation (Figure 3C). The nomogram showed that DNAJB11 levels and pathological stage were independent predictive factors of progression-free interval in THCA patients, while age was not (Figure 3D).

Figure 3.

The prognostic value of the levels of DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11) mRNA in thyroid carcinoma (THCA). (a) Kaplan–Meier survival curve analysis of patients with THCA with respect to the DNAJB11 mRNA levels. The progression-free interval in The Cancer Genome Atlas (TCGA) database cohort was performed for all patients. (b) Low levels of DNAJB11 mRNA were associated with poor progression-free interval in patients with THCA. (c) Calibration plot of the nomogram for predicting the probability of progression-free interval (PFI) at 1, 5 and 10 years. (d) Nomogram for predicting the probability of 1-, 5- and 10-year PFI for patients with THCA. The colour version of this figure is available at: http://imr.sagepub.com.

The Cox regression model was used to analyse the risk factors for the prognosis of patients with THCA. Cox regression analysis showed that DNAJB11 mRNA levels were an independent prognostic factor (HR 0.410; 95% CI: 0.237, 0.766; P = 0.009) in THCA (Table 2). In addition, T stage (HR 2.977; 95% CI 1.032, 8.586; P = 0.044), pathological stage (HR 7.197; 95% CI 2.316, 22.366; P < 0.001) and residual tumour (HR 3.356; 95% CI 1.008, 11.182; P = 0.049) were also independent prognostic factors of PFS in THCA.

Table 2.

Univariate and multivariate Cox regression analysis of risk factors for progression-free survival in patients with thyroid carcinoma.

Characteristic	n	Univariate analysis		Multivariate analysis
Characteristic	n	HR (95% CI)	P-value	HR (95% CI)	P-value
T stage, T3&T4 versus T1&T2	502	2.977 (1.032, 8.586)	P = 0.044	0.288 (0.038, 2.161)	NS
N stage, N1 versus N0	452	1.438 (0.469, 4.407)	NS
M stage, M1 versus M0	291	4.252 (0.908, 19.922)	NS	0.000 (0.000, Inf)	NS
Pathological stage, III&IV versus I&II	502	7.197 (2.316, 22.366)	P < 0.001	11.876 (1.554, 90.748)	P = 0.017
Residual tumour, R1&R2 versus R0	440	3.356 (1.008, 11.182)	P = 0.049	4.226 (0.803, 22.226)	NS
Histological type, classical versus follicular & other & tall cell	502	4.930 (0.649, 37.425)	NS
Age, >45 versus ≤45	502	773760021.346 (0.000, Inf)	NS
Sex, male versus female	502	1.969 (0.712, 5.444)	NS
Race, white versus Asian/black or African-American	410	1.328 (0.301, 5.855)	NS
Neoplasm location, bilateral versus isthmus & left lobe & right lobe	496	1.156 (0.255, 5.235)	NS
Primary neoplasm focus type, multifocal versus unifocal)	492	0.255 (0.058, 1.130)	NS	0.237 (0.028, 2.003)	NS
Thyroid gland disorder history, lymphocytic thyroiditis & nodular hyperplasia & other, specify versus normal)	444	0.182 (0.024, 1.404)	NS
Extrathyroidal extension, yes versus no	484	2.241 (0.831, 6.040)	NS
BRAF status, Mut versus WT	485	0.608 (0.213, 1.738)	NS
NRAS status, Mut versus WT	485	0.000 (0.000, Inf)	NS
HRAS status, Mut versus WT	485	1.911 (0.247, 14.760)	NS
DNAJB11 mRNA levels, high versus low	502	0.410 (0.237, 0.766)	P = 0.009

HR, hazard ratio; CI, confidence interval; T, tumour; N, node; M, metastasis; BRAF, v-raf murine sarcoma viral oncogene homolog B1; MUT, mutant; WT, wild type; NRAS, neuroblastoma RAS viral oncogene homolog; HRAS, Harvey rat sarcoma viral oncogene homolog; DNAJB11, DnaJ Heat Shock Protein Family (Hsp40) Member B11; NS, no significant association (P ≥ 0.05).

To elucidate whether DNAJB11 could play a role in promoting THCA occurrence, RNA sequence gene expression analysis was performed to compare the gene expression profiles of the DNAJB11-high and -low groups. Using the cut-off criteria of adjust P-value <0.05, |logFC| > 1.5, the data from the TCGA database were analysed using the DESeq2 package in R.¹⁷ A total of 177 upregulated mRNA and 88 downregulated mRNAs were detected between the DNAJB11-high and -low groups (Figure 4A). The co-occurrence differential mRNA of the DNAJB11-high and -low groups are illustrated by a Heat map (Figure 4B). To predict the functional enrichment information of interactive mRNA of DNAJB11, clusterProfiler package¹⁹ was used for gene associated DEG in THCA. The top 20 GO enrichment items were classified into three functional groups: biological process group (13 items) (Figure 4C), cellular components (two items) (Figure 4D) and molecular functions analysis (five items) (Figure 4E). The biological processes of DNAJB11 included extracellular structure organization, peptide hormone secretion, extracellular matrix organization, signal release, hormone secretion, hormone transport, cartilage development, multi-multicellular organism process, positive regulation of synapse assembly, female gonad development, regulation of hormone secretion, connective tissue development and development of primary female sexual characteristics. The cellular components for these genes were involved in the collagen-containing extracellular matrix and ER lumen. The molecular functions for DNAJB11 were mainly transcription regulation by receptor ligand activity, cytokine receptor binding, cytokine activity, metal ion transmembrane transporter activity and extracellular matrix structural constituent. The analysis also found that DNAJB11 was involved in immune-related biological process, including neutrophil migration and neutrophil chemotaxis. The DNAJB11 co-occurrence derived protein–protein interaction network was constructed using the Search Tool for the Retrieval of Interacting Genes (Figure 4F).²²

Figure 4.

Identifying genes with differential expression between DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11)-high and -low patients with thyroid carcinoma. (a) Volcano plot of the differentially expressed mRNA. (b) Heat map of the co-occurrence differential mRNA of the DNAJB11-high and -low groups. Bubble graphs showing significantly enriched biological process (c), cellular component (d) and molecular function (e) of gene ontology (GO) annotations of DNAJB11-related mRNA. (f) The protein–protein interaction network of significant DNAJB11 co-occurrence genes. The colour version of this figure is available at: http://imr.sagepub.com.

To identify signalling pathways that were differentially activated in THCA, a GSEA was undertaken between the DNAJB11-low and -high data sets. GSEA revealed significant differences (FDR <0.25, adjust P-value <0.05) in enrichment of Molecular Signature Database Collection (h.all.v7.0.symbols). The results showed that the DNAJB11 low mRNA level phenotype consisted of many key pathways and correlated with tumorigenesis. The analysis showed that coagulation, apoptosis, complement, epithelial mesenchymal transition, mitotic spindle, p53 pathway, androgen response and oestrogen response were differentially enriched in the DNAJB11 low level phenotype.

The association between the level of DNAJB11mRNA and immune cell infiltration level quantified by ssGSEA in the THCA microenvironment was analysed using Spearman's rank correlation coefficient. The results showed that the DNAJB11 mRNA level was significantly inversely correlated with the abundance of innate immunocytes (DCs and T_reg) (Figures 5A and 5B; P < 0.001). Wilcoxon signed rank test also showed that the enrichment scores of DCs, T_reg cells, macrophages, aDCs, iDCs, neutrophils, T_h1, T_h2, mast cells, Tem, cytotoxic cells, T cells, B cells and Tcm cells were significantly higher in DNAJB11-low level samples than in DNAJB11-high level samples (Figure 5C; P < 0.001).

Figure 5.

The correlation between the level of DnaJ Heat Shock Protein Family (Hsp40) Member B11 (DNAJB11) mRNA and immune infiltration in the tumour microenvironment. (a) Correlation between the relative abundances of 24 immune cell types and DNAJB11 mRNA levels. The size of dots shows the absolute value of Spearman's rank correlation coefficient. (b) DNAJB11 mRNA had a significant inverse correlation with dendritic cells (DCs), regulatory T (T_reg) cells and macrophage infiltrating levels. (c) The Wilcoxon signed rank test revealed that the enrichment scores of DCs, T_reg cells and macrophages were significantly higher in low DNAJB11 mRNA samples than in the high DNAJB11 mRNA samples. Data presented as mean ± SD. aDCs, activated dendritic cells; iDCs, immature dendritic cells; T_h1, type 1 helper T cell; T_h2, type 2 helper T cell; Tem, effector memory T cell; Tcm, central memory T cell; NK, natura killer; T_h17, type 17 helper T cell; Tfh, follicular helper T cell; Tgd, gamma delta T; pDCs, plasmacytoid dendritic cell. The colour version of this figure is available at: http://imr.sagepub.com.

Discussion

The ER stress response has been identified as a key driver in the survival of malignant cells and the progression of cancer.¹⁴ DNAJB11 protects against ER stress by stabilizing proteins.¹⁵ The function and mechanisms of DNAJB11 are only now being revealed. For example, DNAJB11 has been implicated in autosomal-dominant polycystic kidney disease.²³ DNAJB11 is known to be closely related to a variety of signalling pathways, including abnormal signalling pathways that are prevalent in tumours.¹⁶ However, little is known about DNAJB11 mRNA levels and their prognostic significance in THCA. To gain a better understanding of the possible functions of DNAJB11 in THCA and its regulatory network, bioinformatics analysis of publicly available sequencing data was used to guide future research in THCA.

According to research, the role of DNAJB11 in tumorigenesis can be either inhibitory or proliferative depending on the specific tissues.²⁴ DNAJB11 promotes the progression of hepatocellular carcinoma by inhibiting alpha-1-antitrypsin mutant Z degradation.¹² The levels of DNAJB11 are increased in oral cavity squamous cell carcinoma tissues when compared with adjacent noncancerous epithelia.¹³ K1 is at the far-left end of the Kaposi sarcoma-associated herpesvirus genome that encodes a unique transmembrane glycoprotein; and DNAJB11 is required for its anti-apoptotic function.²⁵ Overexpression of the GRP78 gene partially suppressed the apoptosis and induction of C/EBP homologous protein in the presence of celecoxib and this suppression was stimulated by co-expression of DNAJB11.²⁶ Based on the TCGA database, this current study found that DNAJB11, compared with normal tissues, was present in significantly higher levels in BLCA, BRCA, CESC, CHOL, COAD, ESCA, GBM, HNSC, KIRC, KIRP, LIHC, LUAD, LUSC, PRAD, READ, SARC, STAD and UCEC, while it was present at significantly lower levels in KICH and THCA. These current findings showed that a lower DNAJB11 level predicted poor THCA outcomes. The prognostic value of all 502 patients with THCA from the TCGA database was evaluated. Previously, the link between DNAJB11 levels and the stage of THCA tissues has not been evaluated. However, in this current study, THCA had lower DNAJB11 levels than normal tissues and this lower DNAJB11 level was associated with an advanced clinical stage. This suggests that DNAJB11 in THCA could act as a tumour suppressor gene. In addition, DNAJB11 level was found to be correlated with poor THCA prognosis in pathological stages I to IV, T2 to T4 and N0 to N1, with the highest HR for poor PFS being associated with a low DNAJB11 level in THCA. Cox proportional hazard analysis also confirmed this result, demonstrating that DNAJB11 is an independent prognostic marker for THCA.

Disruption of ER function leads to abnormal accumulation of unfolded proteins and promotes apoptosis-induced cell death in a process known as ‘ER stress’.²⁷ The ER stress response has been identified as a key driver in the survival of malignant cells and the progression of cancer, and the increased induction of the ER stress pathway has been identified as a hallmark of cancer cells.¹⁴ER stress inhibition has been shown to slow the progression of cancers such as pancreatic carcinoma,²⁸ breast carcinoma²⁹ and multiple myeloma.³⁰ DNAJB11 protects against ER stress by stabilizing proteins.¹⁵ To investigate the function of DNAJB11 in THCA, this current study undertook GSEA using data from the TCGA database. The current results showed that coagulation, apoptosis, complement, epithelial mesenchymal transition, mitotic spindle, p53 pathway, androgen response and oestrogen response are differentially enriched in the DNAJB11 low level phenotype. This suggests that DNAJB11 could be used as a prognostic indicator and therapeutic target in THCA. In addition, this current study evaluated the association between DNAJB11 and immune cell populations. DNAJB11 was found to have a significant correlation with DCs and T_reg cells. DNAJB11 may be required for effective anti-tumour immune responses.

Furthermore, this current study found that DNAJB11 was involved in immune-related biological processes, including neutrophil migration and neutrophil chemotaxis. These analyses showed that the levels of immune infiltration and diverse immune marker sets were correlated with DNAJB11 level in THCA. These current findings demonstrated that there was a moderate-to-strong inverse correlation between DNAJB11 level and DC, T_reg cell and macrophage infiltration. T_reg cells are immunosuppressive immune cells that are involved in the development of tumours.³¹ The presence of T_reg cells in tumour tissues has been shown to impair anti-tumour immunity, implying a poor prognosis for affected patients.³² Suppression of T_reg cell function is considered to be an effective strategy for cancer therapy.³¹ In addition, the immune cells that are effective are suppressed in response to an increase in T_reg cells.³¹ These immune effective cells are crucial in tumour escape.³³ ER stress plays a critical role in the induction of inflammation in a variety of pathological conditions.¹⁴

Recent research indicates that macrophages are primarily responsible for tumour promotion in the immune microenvironment.³⁴ M2-type tumour associated macrophages have been linked to the aggressiveness of various cancers.³⁵ Thus, this current research contributes to a better understanding the potential role of DNAJB11 in THCA immunology and its application as a cancer biomarker.

This current study had several limitations. First, the correlation between DNAJB11 mRNA and protein levels should be validated using cell experiments and clinical samples, as predicting protein levels solely on the basis of mRNA levels is far from perfect. Secondly, larger studies of clinical data should be used to assess the relationship between DNAJB11 levels and clinical features, THCA stage and prognosis.

In conclusion, DNAJB11 may be associated with a poor prognosis and may play important roles in the regulation of the ER stress pathway in THCA. This current study revealed some of the roles of DNAJB11 in THCA and identified a potential biomarker for THCA diagnosis and prognosis.

Footnotes

Author contributions

Rongxin Sun and Longyan Yang were involved in the conception and design of this study. Rongxin Sun, Longyan Yang and Yan Wang performed the data analysis and interpretation of the results. Rongxin Sun, Longyan Yang and Yan Wang prepared the first draft of the manuscript. Yuanyuan Zhang, Jing Ke and Dong Zhao did a critical revision of the manuscript. Dong Zhao supervised the study. All authors read and approved the final manuscript.

Declaration of conflicting interest

The authors declare that there are no conflicts of interest.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and publication of this article: This work was supported by the Luhe Hospital Science and Technology Innovation Achievements Incubation Incentive Special Project (no. LHYY2021-CX02).

ORCID iD

Dong Zhao

References

Siegel

Miller

Jemal

Cancer statistics, 2020. CA Cancer J Clin 2020; 70: 7–30.

Shi

Sun

, et al. Identification and validation of a pyroptosis-related prognostic signature for thyroid cancer. Cancer Cell Int 2021; 21: 523.

Morris

Tuttle

Davies

Changing Trends in the Incidence of Thyroid Cancer in the United States.

JAMA Otolaryngol Head Neck Surg 2016; 142: 709–711.

Fahiminiya

de Kock

Foulkes

WD.

Biologic and Clinical Perspectives on Thyroid Cancer. N Engl J Med 2016; 375: 2306–2307.

Peng

Croce

CM.

The role of MicroRNAs in human cancer. Signal Transduct Target Ther 2016; 1: 15004.

Kuo

Ren

Lao

, et al. Suppression of polyglutamine protein toxicity by co-expression of a heat-shock protein 40 and a heat-shock protein 110. Cell Death Dis 2013; 4: e833.

Shen

Meunier

Hendershot

LM.

Identification and characterization of a novel endoplasmic reticulum (ER) DnaJ homologue, which stimulates ATPase activity of BiP in vitro and is induced by ER stress.

J Biol Chem 2002; 277: 15947–15956.

Jin

Zhuang

Hendershot

LM.

ERdj3, a luminal ER DnaJ homologue, binds directly to unfolded proteins in the mammalian ER: identification of critical residues.

Biochemistry 2009; 48: 41–49.

Shen

Hendershot

LM.

ERdj3, a stress-inducible endoplasmic reticulum DnaJ homologue, serves as a cofactor for BiP's interactions with unfolded substrates.

Mol Biol Cell 2005; 16: 40–50.

10.

Schorr

Klein

Gamayun

, et al. Co-chaperone Specificity in Gating of the Polypeptide Conducting Channel in the Membrane of the Human Endoplasmic Reticulum. J Biol Chem 2015; 290: 18621–18635.

11.

Guo

Snapp

EL.

ERdj3 regulates BiP occupancy in living cells. J Cell Sci 2013; 126: 1429–1439.

12.

Pan

Cao

Gong

The endoplasmic reticulum co-chaperone ERdj3/DNAJB11 promotes hepatocellular carcinoma progression through suppressing AATZ degradation.

Future Oncol 2018; 14: 3001–3013.

13.

Hsu

Peng

, et al. Secretome profiling of primary cells reveals that THBS2 is a salivary biomarker of oral cavity squamous cell carcinoma. J Proteome Res 2014; 13: 4796–4807.

14.

Yadav

Chae

Kim

, et al. Endoplasmic reticulum stress and cancer. J Cancer Prev 2014; 19: 75–88.

15.

Chen

Chowdhury

, et al. The endoplasmic reticulum HSP40 co-chaperone ERdj3/DNAJB11 assembles and functions as a tetramer. EMBO J 2017; 36: 2296–2309.

16.

Guorong

Yunfeng

Yanli

, et al. Expression and clinical significance of DNAJB11 in epithelial ovarian cancer. Zhejiang Da Xue Xue Bao Yi Xue Ban 2017; 46: 173–178 [Article in Chinese, English abstract].

17.

Love

Huber

Anders

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.

Genome Biol 2014; 15: 550.

18.

Zhou

Pache

, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun 2019; 10: 1523.

19.

Wang

Han

, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012; 16: 284–287.

20.

Bindea

Mlecnik

Tosolini

, et al . Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 2013; 39: 782–795.

21.

Liu

Lichtenberg

Hoadley

, et al. An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics. Cell 2018; 173: 400–416.e11.

22.

Szklarczyk

Gable

Lyon

, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 2019; 47: D607–D613.

23.

Cornec-Le Gall

Olson

Besse

, et al. Monoallelic Mutations to DNAJB11 Cause Atypical Autosomal-Dominant Polycystic Kidney Disease. Am J Hum Genet 2018; 102: 832–844.

24.

Liu

Rios

, et al. Heat shock proteins and cancer. Trends Pharmacol Sci 2017; 38: 226–256.

25.

Wen

Damania

Hsp90 and Hsp40/Erdj3 are required for the expression and anti-apoptotic function of KSHV K1.

Oncogene 2010; 29: 3532–3544.

26.

Tsutsumi

Namba

Tanaka

, et al. Celecoxib upregulates endoplasmic reticulum chaperones that inhibit celecoxib-induced apoptosis in human gastric cells. Oncogene 2006; 25: 1018–1029.

27.

Lai

Teodoro

Volchuk

Endoplasmic reticulum stress: signaling the unfolded protein response. Physiology 2007; 22: 193–201.

28.

Gupta

McGrath

Cavener

DR.

PERK regulates the proliferation and development of insulin-secreting beta-cell tumors in the endocrine pancreas of mice.

PLoS One 2009; 4: e8008.

29.

Nagelkerke

Bussink

Mujcic

, et al. Hypoxia stimulates migration of breast cancer cells via the PERK/ATF4/LAMP3-arm of the unfolded protein response. Breast Cancer Res 2013; 15: R2.

30.

Pike

Singleton

Buffa

, et al. Transcriptional up-regulation of ULK1 by ATF4 contributes to cancer cell survival. Biochem J 2013; 449: 389–400.

31.

Hossain

Liu

Dai

, et al. Reinvigorating exhausted CD8+ cytotoxic T lymphocytes in the tumor microenvironment and current strategies in cancer immunotherapy. Med Res Rev 2021; 41: 156–201.

32.

Deng

Lira

Hudson

, et al. Recombinant Listeria promotes tumor rejection by CD8+ T cell-dependent remodeling of the tumor microenvironment. Proc Natl Acad Sci U S A 2018; 115: 8179–8184.

33.

Mougiakakos

Johansson

Trocme

, et al. Intratumoral forkhead box P3-positive regulatory T cells predict poor survival in cyclooxygenase-2-positive uveal melanoma. Cancer 2010; 116: 2224–2233.

34.

Duan

Luo

Targeting macrophages in cancer immunotherapy. Signal Transduct Target Ther 2021; 6: 127.

35.

Ichimura

Abe

Morikawa

, et al. Low density of CD204-positive M2-type tumor-associated macrophages in Epstein-Barr virus-associated gastric cancer: a clinicopathologic study with digital image analysis. Hum Pathol 2016; 56: 74–80.