Clinical prognosis and bioinformatic analysis of primary thyroid lymphoma.

ABSTRACT: Primary thyroid lymphoma (PTL) is a rare malignant disease with the most common histological type of diffuse large B-cell lymphoma (DLBCL). Hashimoto's thyroiditis (HT) is closely related to the pathogenesis of PTL. The present study is to explore the clinical prognosis of PTL and analyze the gene correlations between PTL and HT.Thirty-nine patients diagnosed with PTL between 2010 and 2018 in our institute were retrospectively reviewed and clinical features were evaluated on PTL survival. Then, overlapping differentially expressed genes (DEGs) between PTL and HT were evaluated for gene ontology, pathways enrichment, protein-protein interaction network analysis. Furthermore, we used gene expression profiling interactive analysis to evaluate the differential expression of these hub genes.In this analysis, International Prognostic Index (IPI) score ≥3 and high β2-MG (>3 mg/L) were associated with worse prognosis in PTL. Notably, a total of 15 both upregulated DEGs in DLBCL and HT were identified and 10 hub genes with a high degree of connectivity were picked out. Among these 10 hub genes, IL6, IL10, CXCL10, and CXCR3 were higher expressed in DLBCL than the normal tissue but have no significant prognosis of DLBCL.High IPI score and high β2-MG level have a poor prognosis in PTL. Besides, IL6, IL10, CXCL10, and CXCR3 are associated with both DLBCL and HT and may be used for the early diagnosis of PTL.

Introduction

Primary thyroid lymphoma (PTL) is a very rare thyroid malignancy that accounts for only 1% to 5% of all thyroid malignancies.[ The most common histopathological type is diffuse large B-cell lymphoma (DLBCL), which accounts for 50% to 70% of cases. Hashimoto's thyroiditis (HT) is an autoimmune condition that leads to the destruction of thyroid cells and is the leading cause of hypothyroidism worldwide.[ A growing amount of evidence indicates that PTL is associated with a background of HT.[ Patients with HT are 40 to 80 times more at risk to develop lymphoma; however, only 0.6% of those with thyroiditis develop lymphoma.[ Moreover, some authors even suggest that all PTLs originate in an HT setting.[ PTL is presumed to develop from lymphoid tissue acquired during the course of an autoimmune process. Nearly 100% PTL patients have HT as revealed by thorough histopathological examination.[ Although PTL is closely related to HT in the process of its occurrence and development, the mechanisms underlying this association are uncertain. One possible mechanism may be that chronic inflammation caused by oxidative stress is involved in the pathogenesis of HT and PTL. Oxidative stress is defined as a cellular imbalance of oxidants over antioxidants.[ Continuous oxidative stress can lead to chronic inflammation, which in turn can mediate most chronic diseases, including lymphoma.[ Studies have found that chronic inflammation is associated with the risk of non-Hodgkin's lymphoma and that the imbalance between reactive oxygen species and antioxidant defense mechanisms leads to increased oxidative stress in autoimmune thyroiditis.[ Therefore, finding new diagnostic targets or oxidative stress markers is of great significance for the early diagnosis of PTL cancer.

To identify genetic changes associated with cancer, bioinformatics analysis has been commonly applied in the field of cancer research. Differentially expressed genes (DEGs) in various types of cancer have been identified by bioinformatics analysis so as to determine their roles in biological processes, molecular functions, and different pathways.[ Therefore, we screened four public mRNA datasets through the bioinformatics tools to identify the correlated genes between DLBCL and HT. Subsequently, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses and a protein-protein interaction (PPI) network analysis were performed to assist in demonstrating the molecular pathogenesis underlying the carcinogenesis and development of PTL.

Moreover, we further explore the clinical prognosis of PTL and analyze the common biomolecular profiles between HT and PTL and evaluate the effect of these co-expressed genes on prognosis.

Materials and methods

Patient characteristics and ethics statement

Thirty-nine patients treated for PTL between 2010 and 2018 at our institution were retrospectively analyzed. Collection of blood for analysis occurred at diagnosis, receipt of any treatment. Clinical data were reviewed, including age, gender, stage, histopathologic type, the presence of HT, diagnostic methods, and treatment types. All patients in the study were followed up, and their recurrence, metastasis, and death were recorded. The follow-up period continued through December 31, 2018. The study was exempted from informed consent and ethical approval was obtained from the ethical committee of our hospital (Ethics number: 2020-L051).

Microarray data and DEGs identification

GSE74266 and GSE29315 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). In total, 29 activated B-cell-like (ABC) DLBCL samples and 29 germinal center B-cell like (GCB) DLBCL samples from GSE74266 and 6 HT samples from GSE29315 and 8 normal samples were selected. Then, GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r), a program that allows the investigator to compare two or more groups of samples within a GEO dataset, was applied to screen the DEGs between DLBCL and HT.[ The log-fold change (FC) in expression and adjusted P-values (adj. P) were determined. The adj. P using the Benjamini-Hochberg method with default values were applied to correct the potential false-positive results. Genes that met the specific cut-off criteria of adj. P < 0.05 and |logFC|>2 were regarded as DEGs. The Venn diagram was performed by using R VennDiagram package.

GO annotation and KEGG pathway enrichment analyses of DEGs

The functions and pathway enrichment of candidate DEGs were analyzed using the WebGestalt database (http://www.webgestalt.org/option.php). Then, enrichment data from the WebGestalt database were extracted using the R Bioconductor ggplot2 and dplyr packages (http://bioconductor.org/biocLite.R) and the bubble plot was applied to visualize the enrichment result.

PPI network construction and hub gene identification

To obtain a PPI map, Search Tool for the Retrieval of Interacting Genes (STRING; http://STRING-db.org) (version 11.0) online database[ was applied to construct the PPI network of 15 co-DEGs. Then, we utilized Cytoscape (v3.7.2) plugin cytoHubba[ to visualize the protein interaction relationship network and analyzed hub genes, which were important nodes with many interaction partners. In this experiment, the genes with the top 10 highest degree values were considered as hub genes.

Hub gene survival analysis

The validation of the expression of the core hub genes in lymphoma and normal tissues was performed and visualized in a boxplot using the online tool gene expression profiling interactive analysis (GEPIA). Then, GEPIA survival assessment was used to investigate the overall association with survival of 4 hub genes from both the low and high-expression lymphoma groups.

Statistical analysis

The survival curves were plotted according to the Kaplan-Meier method and the log-rank test was performed to test the statistical significance. Statistical analyses were performed using SPSS 25.0 (IBM Corp, Armonk, NY, USA). P < 0.05 indicated significant differences. Progression-free survival (PFS) was defined as the time from the treatment until the progression of lymphoma or death. Overall survival (OS) was defined as the time span from the date of diagnosis to the date of death due to any cause.

Results

Patients’ characteristics

The baseline clinical parameters of the patients are presented in Table 1. The patients’ median age was 59 years (range 31–82 years old). The median follow-up time was 49 months (7–134 months). The median treatment cycle was 6 cycles (4–8 cycles). All the PTL patients were diagnosed with DLBCL and none with mucosa-associated lymphoid tissue lymphoma. Using the Lugano staging system, 11 patients (28.2%) were diagnosed with Stage III and IV. According to the algorithm of Hans et al,[ 16 patients (53.3%) were assigned as non-GCB phenotype. Immunohistochemical staining and scoring for CD10, Bcl-6, MUM-1, Bcl-2, and Ki67 were performed on formalin-fixed paraffin-embedded tissues.

Table 1

The clinical characteristics of the 39 patients of PTL.

Characteristics	Number of cases (%)
Age (years)
>60	17 (43.6%)
Female	32 (82.1%)
IPI score of 3-5	5 (12.8%)
Stage III-IV	11 (28.2%)
COO (Hans)
Non-GCB	16 (53.3%)
TgAb positive	13 (40.6%)
TPOAb positive	18 (56.3%)
Immunohistochemical features
Ki67 expression >70%	24 (66.7%)
Bcl-6 positive expression	20 (76.9%)
CD10 positive expression	10 (37.0%)
MUM1 positive expression	14 (46.7%)
Bcl-2 positive expression	9 (60.0%)
ESR (>20 mm/h)	18 (66.7%)
Serum LDH (>245 U/L)	25 (64.1%)
β2-MG (>2.2 mg/L)	13 (33.3%)
Hypothyroidism	10 (31.3%)

β2-MG = β2-microglobulin, COO = cell of origin, ESR = erythrocyte sedimentation rate, GCB = germinal center B-cell, IPI = International Prognostic Index, LDH = lactate dehydrogenase, TgAb = thyroglobulin antibody, TPOAb = thyroid peroxidase antibody.

Prognostic factors analysis for PTL patients

The median OS and PFS time of patients with PTL was 51 months (range, 9–138 months) and 43 months (range, 1–138 months), respectively. A univariate analysis identified several prognostic factors for survival in patients with PTL. High β2-GM level (>2.2 mg/L) patients have decreased OS compared with the normal groups (P = 0.021) (Fig. 1E). Moreover, IPI ≥3 scores significantly affected survival (P < 0.001) (Fig. 1B). However, the survival rate was not significantly influenced by COO, ESR (>20 mm/h), serum lactate dehydrogenase (LDH) (>245 U/L), and HT status (Fig. 1A, C, D, and F). To extend our findings, IPI ≥3 scores (Fig. 1H) and β2-MG (>2.2 mg/L) (Fig. 1K) affected PFS respectively, while COO, ESR (>20 mm/h), serum LDH (>245 U/L), and HT status showed no significant differences in PFS (Fig. 1G, I, J, and L).

Figure 1

Prognostic indicators for OS and PFS in PTL patients. Univariate analysis demonstrated that (B) IPI score 3, (E) high β2-MG level were risk factors associated with OS. Moreover, (H) IPI score ≥3, (K) high β2-MG level were risk factors associated with PFS in PTL patients. β2-MG = β2-microglobulin, COO = cell of origin, ESR = erythrocyte sedimentation rate, GCB = germinal center B-cell, IPI = International prognostic index, LDH = lactate dehydrogenase, OS = overall survival, PFS = progression-free survival.

Immunohistochemical survival analysis

In terms of immunohistochemistry (Fig. 2), only 37.0% (10/27) expressed CD10, 76.9% (20/26) Bcl-6, 46.7% (14/30) MUM1, and 60% (9/15) Bcl-2 and the expression of proliferation index Ki67 ranged from 30% to 90%. According to the Hans classification, 14 patients (46.7%) were germinal center B-cell (GCB) type, 16 patients (53.3%) were non-GCB type. Then, the expression of each immunophenotype biomarker and its relationship with survival were obtained by performing univariate analysis. Our result showed that CD10, Bcl-6, MUM1, Bcl-2, or Ki67 were not statistically significant with OS.

Figure 2

Immunohistochemical survival analysis of PTL patients. (A) Histopathology showed a number of single lymphoid cells infiltrated the thyroid tissue (HE staining ×200). Immunohistochemical staining of PTL. (B) Negative CD10 protein expression (×200). (C) Positive Bcl-6 protein expression (×200). (D) Positive MUM1 protein expression (×200). (E) Negative Bcl-2 protein expression (×200). (F) Ki-67 diffuse, intense positivity (×200).

Identification of DEGs

In this study, we selected two groups of gene expression profiles (GSE74266 and GSE29315). Via GEO2R, 202 DEGs (105 upregulated and 97 downregulated), and 470 DEGs (301 upregulated and 169 downregulated) were identified from the expression profile database GSE74266 and GSE29315, respectively. Then, DLBCL, DEGs, and HT (relative to the normal group) DEGs were subjected to Venn analysis to obtain Venn analysis result (Fig. 3A).

Figure 3

(A) Venn diagram shows 15 DEGs that are both upregulated between G1 up and G2 up. The bubble diagrams display the enrichment results of the 15 upregulated DEGs between DLBCL and HT. (B) Biological processes. (C) Molecular functions. (D) Cellular components. (E) KEGG pathway analysis. (F) Protein-protein interaction network was constructed with 23 nodes and 106 edges using STRING database. (G) The top 10 hub genes such as IL10RA, IL6, STAT3, IL10, IL4, CXCL10, CXCR3, CCL5, GZMB, and PRF1 were screened by Cytoscape (v3.7.2) plugin cytoHubba. DEGs = differentially expressed genes, G1 down = downregulated DEGs in DLBCL, G1 up = upregulated DEGs in DLBCL, G2 down = downregulated DEGs in HT, G2 up = upregulated DEGs in HT.

GO annotation and KEGG pathway enrichment analyses

To obtain a deeper insight into the biological roles of these 15 both upregulated DEGs, the WebGestalt database (http://www.webgestalt.org/option.php) was employed to conduct GO annotation and KEGG pathway enrichment analyses. Then we used the R Bioconductor ggplot2 and dplyr packages (http://bioconductor.org/biocLite.R) to visualize the enrichment result using the bubble plot. Figure 3B-E lists the most relevant enriched GO terms and KEGG pathways. GO BP analysis revealed that the 15 upregulated DEGs were markedly enriched in the immune effector process, the response to other organisms, the external biotic stimulus-response and biotic stimulus-response (Fig. 3B). For GO MF analysis, the top significantly enriched term was the non-membrane spanning protein tyrosine kinase activity (Fig. 3C). The top four significantly enriched CC terms included lysosome, lytic vacuole, nuclear nucleosome, and the integral component of the nuclear inner membrane (Fig. 3D). In addition, the top four markedly enriched pathways for these both upregulated DEGs were graft-versus-host disease, transcriptional misregulation in cancer, cytosolic DNA-sensing pathway, and acute myeloid leukemia (Fig. 3E).

PPI network of DEGs and hub gene identification

In this study, the STRING tool was applied to identify the PPI networks of DEGs. Then, we utilized Cytoscape to visualize the PPI network that consisted of 23 nodes and 106 edges (Fig. 3F). Moreover, the top 10 degree genes in the PPI network were selected as hub genes using cytoHubba, eg IL10RA, IL6, STAT3, IL10, IL4, CXCL10, CXCR3, CCL5, GZMB, and PRF1 genes (Fig. 3G).

CXCL10, CXCR3, IL6, and IL10 might be the diagnostic biomarkers for PTL

Moreover, we performed the differential expression analysis using GEPIA tools to further clarify the expression of these hub genes in DLBCL and normal tissues. From the box-profiles in Fig. 4A–D, these results indicated that 4 genes in the 10 hub genes were expressed at a higher level in DLBCL than the normal tissues (p < 0.05), for example, CXCL10, CXCR3, IL6, and IL10 genes. In addition, we detected the relationship between the above 4 hub genes and the survival of DLBCL patients. As is shown in Fig. 4E-H, CXCL10, CXCR3, IL6, and IL10 hub genes have no significant prognosis of DLBCL (P < 0.05). Our results suggest that CXCL10, CXCR3, IL6, and IL10 genes can be used for the early diagnosis of PTL.

Figure 4

Validation of the mRNA expression levels of (A) CXCL10, (B) CXCR3, (C) IL6, and (D) IL10 in DLBC tissues and normal tissues using GEPIA. These four box plots are based on 47 DLBC samples (marked in red) and 337 normal samples (marked in gray). ∗P < 0.05 was considered statistically significant. DLBC, diffuse large B-cell lymphoma. The four hub genes (E) CXCL10, (F) CXCR3, (G) IL6, and (H) IL10 impacted the survival rate of DLBCL patients. Log-rank P < 0.01 was regarded as statistically significant.

Discussion

PTL, a rare thyroid malignant tumor, is defined as a lymphoma involving only the thyroid and or regional lymph nodes without metastasis to other areas upon diagnosis.[ The most common histological type of PTL is DLBCL, which is considered a high-grade lymphoma with a more aggressive clinical course than other subtypes of lymphoma.[ Moreover, HT is the most important risk factor and it causes a 40- to 80-fold increase in the risk of PTL.[ However, the process that HT developed into PTL is not clear. Although a combination of chemotherapy and locoregional radiation is recognized as an effective treatment for PTL, 5-year OS purview of PTL patients is still between 50% and 60%.[ Thus, identifying novel diagnostic and therapeutic markers for PTL continues to be a significant research focus. Bioinformatics analysis has been widely used to investigate genetic alterations in the progression of diseases and may enable us to identify novel diagnostic and therapeutic targets.[ Since DLBCL accounts for the most common histological type of PTL, we screened common biomarkers between HT and DLBCL by applying bioinformatics analysis and identified the potential prognostic factors for PTL. In the present study, the medical records of 39 PTL patients in our institution were retrospectively reviewed and survival analysis indicated that high IPI scores and high β2-MG level were associated with worse prognosis in PTL. Then 202 DLBCL DEGs (105 upregulated and 97 downregulated), and 470 HT DEGs (301 upregulated and 169 downregulated) were identified from the expression profile database GSE74266 and GSE29315, respectively. At last, a total of 15 both upregulated DEGs between DLBCL and HT were revealed by applying GEO2R tools. To further understand these DEGs, we performed GO function and KEGG pathways analyses. GO BP annotation indicated that the immune effector process, the response to other organisms, the external biotic stimulus-response and biotic stimulus-response were markedly enriched. Besides, the top significantly enriched term was the non-membrane spanning protein tyrosine kinase activity in GO MF analysis, and the top four significantly enriched CC terms included lysosome, lytic vacuole, nuclear nucleosome, and integral component of the nuclear inner membrane. Moreover, KEGG pathway showed that the upregulated DEGs included graft-versus-host disease, transcriptional misregulation in cancer, cytosolic DNA-sensing pathway, and acute myeloid leukemia. These results revealed that the upregulated DEGs are mainly related to the immune effector process between HT and thyroid lymphoma. And the previous study has shown that thyroid cancer related immune cells can exert both antitumor and protumor functions in thyroid cancer.[

To further screen hub genes, a PPI network consisting of 23 nodes and 106 edges was constructed, and 10 hub genes were identified, including IL10RA, IL6, STAT3, IL10, IL4, CXCL10, CXCR3, CCL5, GZMB, and PRF1 genes. Subsequently, the GEPIA webserver was used to assess the association between hub genes expression and DLBCL prognosis. IL6, IL10, CXCL10, and CXCR3 genes were more highly expressed in DLBCL than the normal tissues but the OS analysis indicated that high expression of these genes was not significantly related to the poor survival of DLBCL patients.

IL6 is an inflammation cytokine, which plays an important role in acquired immune response by stimulation of antibody production. Earlier studies had regarded it as a malignancy predictor with relatively high sensitivity and specificity.[ Moreover, in a chronic inflammatory process, IL-6 may generate free radicals that can induce DNA damages, causing tumor initiation.[ IL10, a highly pleiotropic cytokine, is produced by a variety of cells. It also participated in the pathogenesis and development of autoimmune diseases and tumorigenesis.[ Additionally, previous literature proposed that the IL-6 or IL-10 cytokine not only induced oxidative stress but was also associated with a tumorigenic process or poor prognosis.[ Moreover, earlier studies had considered that oxidative stress not only played a crucial role in the biology of lymphomas but was also implicated in the pathogenesis of autoimmune thyroiditis.[ Song et al[ verified that oxidative stress interacted with nuclear transcription factors to promote the expression of inflammatory cytokines. In this study, we found that IL6 was both upregulated in HT and DLBCL samples, and compared with the normal tissues, the level of IL6 significantly increased in DLBCL. Additionally, the level of IL10 significantly increased in the DLBCL group but was not associated with OS in DLBCL. Therefore, we speculate that oxidative stress may be involved in the pathogenesis of HT and PTL via the inflammatory cytokines, which are produced by chronic inflammation.

Previous study indicated that CXCR3 and its ligand CXCL10 influence the tumor microenvironments by regulating angiogenesis, recruiting activated immune cells, and affecting tumor cells by promoting or inhibiting tumor progression.[ Besides, the activation of CXCL10/CXCR3 may produce an inflammatory reaction by generating reactive oxygen species (ROS).[ When ROS was in an excessive production, it may cause oxidative stress. Based on bioinformatics analysis, we found that the level of CXCR3 and its ligand CXCL10 were significantly highly expressed in lymphoma tissue. Our finding agreed with previous studies and it may further indicate that oxidative stress is involved in the formation of PTL through chronic inflammation.

In conclusion, although the clinical data analysis revealed that IPI score ≥3 and high β2-MG (>3 mg/L) were associated with worse prognosis in PTL, our systematic bioinformatics assessment demonstrated that IL6, IL10, CXCL10, and CXCR3 genes might serve as new diagnostic targets during the process of HT to PTL. Besides, oxidative stress may participate in the process of lymphomagenesis via chronic inflammation. Nevertheless, there still exist limitations of our study, since the hub genes we identified are only based on the results of bioinformatics analysis and further experimentations are still needed to determine trusted cut-off values of IL6, IL10, CXCL10, and CXCR3 so as to establish a direct relationship between HT and PTL.

Author contributions

Conceptualization: Zhimin Bai, Lingyu Li.

Data curation: Jiangtao Wang, Jin Zhao.

Formal analysis: Zhimin Bai, Tao Guan.

Investigation: Jiangtao Wang, Jin Zhao.

Methodology: Tao Guan.

Supervision: Liping Su.

Writing – original draft: Zhimin Bai.

Writing – review & editing: Liping Su.