Research Paper Volume 14, Issue 11 pp 4739—4754

A novel immune-related microRNA signature for prognosis of thymoma

Bin Wang1,2, , He Xiao2, , Xin Yang3, , Ying Zeng4, , Zhimin Zhang5, , Rui Yang1,2, , Hang Chen1,2, , Chuan Chen2, *, , Junxia Chen1, *, ,

  • 1 Department of Cell Biology and Genetics, Chongqing Medical University, Chongqing 400016, China
  • 2 Department of Oncology, Daping Hospital, Army Medical University, Chongqing 400042, China
  • 3 Department of Pathology, Daping Hospital, Army Medical University, Chongqing 400042, China
  • 4 Department of Pathology, Chongqing General Hospital, University of Chinese Academy of Sciences, Chongqing 410309, China
  • 5 Department of Thoracic Surgery, The Third Medical Center of PLA General Hospital, Chongqing 100039, China
* Equal contribution

Received: March 6, 2022       Accepted: May 10, 2022       Published: June 7, 2022
How to Cite

Copyright: © 2022 Wang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Introduction: Immune microenvironment and microRNAs serve as common predictors for diagnosis and prognosis of tumors.

Methods: Expression of 122 genes and 126 microRNAs in thymoma was obtained from TCGA database. The proportion of tumor-infiltrating cells was calculated, and IMRS was constructed. TREM2hi score was calculated before functional enrichment analysis on gene sets.

Results: IMRS3, TREM2hi score, and CD8+ T lymphocyte abundance were significantly different among WHO classifications. WHO classification, Masaoka staging, and miR-130b-5p, miR-1307-3p, miR-425-5p, CD8, CD68, and CCL18 expression were prognostic factors for relapse-free survival and overall survival. IMRS3 upregulation polarized macrophages into M2, which rejected CD8+ T and other effector lymphocytes to promote thymoma malignant progression.

Conclusions: BRRS may present a novel immune-related microRNA signature for TET prognosis.


Thymic epithelial tumors (TETs) originate from thymic epithelial cells with a relatively low incidence of approximately 0.13/100,000 [1]. Currently, the recurrence and survival prognosis of TETs mainly rely on the World Health Organization (WHO) classification and Masaoka staging [2]. According to 2015 WHO classification, TETs can be divided into thymoma (types A, AB, B1, B2, and B3) and thymic carcinoma (type C) [3]. Types A, AB, B1, B2, and B3 TETs present a low malignancy, less recurrence, and good survival prognosis. In contrast, type C TETs are considered highly malignant and easy to relapse, together with a poor survival prognosis. In Masaoka staging, stages I-II are deemed to have a low risk of recurrence and a good survival prognosis, while stages III-IV present a high risk of recurrence and a poor survival prognosis. However, there are certain limitations to predicting recurrence and survival prognosis of TET patients as per the WHO classification and Masaoka staging system in a clinical setting. For example, recurrence has also been noted for type A, AB, or B1 thymomas with a low malignancy; early metastasis, among other factors, affects the prognosis of patients [4]. As demonstrated in recent studies, biomarkers that are related to the immune microenvironment are associated with the prognosis of TETs. Therefore, the identification of more accurate biomarkers carries tremendous clinical implications for predicting the prognosis of TETs [5].

Significant differences have been found in the expression of various microRNAs related to the immune microenvironment between normal and tumor cells. These differentially expressed microRNAs can act as proto-oncogenes or tumor suppressor genes by regulating different target genes, and their roles are closely related to the occurrence, development, treatment, and prognosis of tumors [6]. Tumor-associated macrophages are the most abundant immune cells in the tumor microenvironment. MicroRNAs regulate macrophage differentiation, functional polarization, and cellular crosstalk [7]. Abnormally expressed microRNAs can regulate the differentiation, proliferation, apoptosis, and function of immune cells, thereby affecting the occurrence and development of tumors [8]. Moreover, CCL18 is a marker of M2 macrophage subpopulations, and its elevated level in serum and tumors is associated with a poor prognosis [9]. The tumor microenvironment is an essential factor for tumor immune escape. The inhibition of immune response by the tumor microenvironment remains a major clinical concern during the administration of immunotherapy. The regulatory mechanism of microRNAs on immune cells in the tumor microenvironment can provide new targets for tumor immunotherapy [10]. However, few studies have been done on the impact of microRNAs on the immune microenvironment and the prognosis of TET patients. Here, we developed BRRS for TETs using microRNA-mRNA pairs and tumor-infiltrating immune cells from the TCGA dataset and performed validation in an independent cohort.


Identification of microRNAs influencing macrophage function on progression-free interval

In the TCGA thymoma datasets, PFI (progression-free interval) data and cell abundances inferred from EPIC (electronic Publication Information Center) were available for 97 tumor samples. Macrophages were the most pronounced cell population for the prognosis of PFI (Supplementary Table 2). Based on the macrophage abundance obtained by EPIC, 16 microRNAs and 850 genes were identified as significantly interacting with macrophages to affect PFI (Supplementary Table 3). A total of 41 microRNA-mRNA pairs were retrieved from the anamiR analysis (Supplementary Table 4). From the expression levels of these microRNAs supported by the microRNA datasets, LASSO (least absolute shrinkage and selection operator) returned a model containing four microRNAs with non-zero coefficients (Table 1). The baseline characteristics of the 119 samples with PFI data are shown in Supplementary Table 5 (sheet ‘Baseline of TCGA TYHM’). Although WHO histology (C vs. A-AB), hsa-miR-1307-3p, hsa-miR-425-5p, and TREM2hi signature were significantly associated with PFI, the effects of IMRS3 and IMRS4 were more conspicuous (Supplementary Table 5, sheet ‘Univariate Cox regression’). Moreover, IMRS3 and IMRS4 were independent prognostic factors for PFI after adjustment for gender, age, WHO classification, Masaoka stage, and other individual biomarkers (Table 2). Taking HR values and 95% CIs into consideration, IMRS3 may be as reliable as IMRS4 in terms of the prognosis of PFI.

Table 1. MicroRNAs included into model resulting from LASSO.


Table 2. Multivariate Cox regression for PFI illustrating independent prognosis of IMRS3 and IMRS4.

HR (95% CI)PHR (95% CI)P
Age1.004 (0.969-1.041).8271.005 (0.969-1.042).791
Gender (Male vs Female)0.737 (0.288-1.888).5250.720 (0.282-1.835).491
Myasthenia gravis (Yes vs No)0.671 (0.225-2.003).4740.660 (0.221-1.970).456
WHO histology (B1-B3 vs A-AB)1.675 (0.503-5.581).4011.789 (0.538-5.946).343
WHO histology (C vs A-AB)1.039 (0.174-6.222).9660.981 (0.159-6.050).983
Masaoka's stage (III-IVB vs I-II)0.924 (0.265-3.225).9010.858 (0.245-3.003).811
TREM2hi1.334 (0.733-2.429).3461.310 (0.714-2.405).383
IMRS4/IMRS32.633 (1.226-5.657).0132.921 (1.273-6.701).011

Characteristics of immune microenvironment in thymoma

As shown in Figure 1, there was a clear trend that IMRS3 and the TREM2hi signature increased with the increase in the WHO grades. However, the abundance of CD8 T cells exhibited a parabolic change. Detailed pairwise comparisons are shown in Supplementary Table 6. Moreover, there was a significant positive correlation between IMRS3 and the TREM2hi signature (Figure 2). In addition, CCL18 (Chemokine (C-C motif) ligand 18), one of the components of the TREM2hi signature, was also highly positively correlated with miR-1307-3p, miR-425-5p, IMRS3, and TREM2hi signature itself (Figure 2). Gene set enrichment analysis showed that tumors with high TREM2hi signature had several critical immune pathways in addition to “interleukin-10 signaling”, such as “interferon-gamma signaling”, “antigen processing-cross presentation”, and “interferon-alpha/beta signaling” (Supplementary Table 7). These results indicated a complex regulatory network of infiltrating monocyte differentiation in thymoma tissues. Importantly, IMRS3 as well as macrophages, especially TREM2hi macrophages, may be essential for thymoma progression.

Distribution of IMRS3, TREM2hi signature, and abundance of CD8+ T cells among the WHO classification. Data were obtained from 122 thymoma gene expression values, as well as phenotype files and survival data all from the UCSC, expression matrix in the format of log2(x+1) transformed RSEM normalized counts.

Figure 1. Distribution of IMRS3, TREM2hi signature, and abundance of CD8+ T cells among the WHO classification. Data were obtained from 122 thymoma gene expression values, as well as phenotype files and survival data all from the UCSC, expression matrix in the format of log2(x+1) transformed RSEM normalized counts.

Pairwise correlation of TREMhi signature, microRNAs, and eleven genes in TREM2hi signature of 120 THYM samples.

Figure 2. Pairwise correlation of TREMhi signature, microRNAs, and eleven genes in TREM2hi signature of 120 THYM samples.

Correlation of microRNAs and immune-related biomarkers with clinical features

To further explore the roles of miR-130b-5p, miR-1307-3p, miR-425-5p, and the TREM2hi signature in determining RFS of thymoma, these biomarkers were evaluated in an independent cohort of 99 patients. Baseline characteristics of the cohort are shown in Supplementary Table 8. Consistent with the TCGA thymoma dataset, a significantly higher proportion of patients with WHO type C or Masaoka stage III-IV thymoma had increased expression of CD68 and CCL18 (Table 3). Representative IHC images of CD8, CD68, and CCL18 are shown in Figure 3. Furthermore, the expression levels of miR-130b-5p, miR-1307-3p, and miR-425-5p were all markedly increased in WHO type C compared with A-AB or B1-B3 types (Table 4). These results confirmed that miR-130b-5p, miR-1307-3p, and miR-425-5p might be involved in the progression of thymoma.

Table 3. Comparison of expression of CD8, CD68 and CCL18 in subtypes of WHO classification.

Female12 (24.5)37 (75.5)1.081/0.29933 (67.3)16 (32.7)1.846/0.17429 (59.2)20 (40.8)0.842/0.359
Male17 (34.0)33 (66.0)27 (54.0)23 (46.0)25 (50.0)25 (50.0)
Myasthenia gravis
No15 (28.8)37 (71.2)0.011/0.91830 (57.7)22 (42.3)0.389/0.53328 (53.8)24 (46.2)0.022/0.883
Yes14 (29.8)33 (70.2)30 (63.8)17 (36.2)26 (55.3)21 (44.7)
WHO classification
A-AB1 (4.5)21 (95.5)34.372/<0.00120 (90.9)2 (9.1)3.270/<0.00117 (77.3)5 (22.7)21.541/<0.001
B1-B313 (21.7)47 (78.3)39 (65.0)21 (35.0)36 (60.0)24 (40.0)
C15 (88.2)2 (11.8)1 (5.9)16 (94.1)1 (5.9)16 (94.1)
Masaoka's stage
I-II1 (1.7)59 (98.3)56.123/<0.00154 (90.0)6 (10.0)55.118/<0.00149 (81.7)11 (18.3)45.186/<0.001
III-IV28 (71.8)11 (28.2)6 (15.4)33 (84.6)5 (12.8)34 (87.2)
Adjuvant radiotherapy
No26 (32.5)54 (67.5)2.070/<0.15047 (58.8)33 (41.3)0.601/0.43842 (52.5)38 (47.5)0.703/0.402
Yes3 (15.8)16 (84.2)13 (68.4)6 (31.6)12 (63.2)7 (36.8)
Adjuvant chemotherapy
No28 (33.7)55 (66.3)3.655/0.05648 (57.8)35 (42.2)1.656/0.19842 (50.6)41 (49.4)3.220/0.073
Yes1 (6.3)15 (93.8)12 (75.0)4 (25.0)12 (75.0)4 (25.0)
Representative IHC images showing different expression of CD8 (A1–A3), CD68 (B1–B3), and CCL18 (C1–C3). A1, B1, and C1 indicates + IHC score; A2, B2, and C2 indicates ++ IHC score; A3, B3, and C3 indicates +++ IHC score. (n = 3 for A–C groups; scale bars, 200 μm).

Figure 3. Representative IHC images showing different expression of CD8 (A1A3), CD68 (B1B3), and CCL18 (C1C3). A1, B1, and C1 indicates + IHC score; A2, B2, and C2 indicates ++ IHC score; A3, B3, and C3 indicates +++ IHC score. (n = 3 for AC groups; scale bars, 200 μm).

Table 4. Comparison of expression level of microRNAs among WHO classification, Masaoka’s stage and CD8, CD68 and CCL18 lowly and highly expressed subgroups determined by immunohistochemistry.

WHO classification
A-AB0.085 (-0.610-0.261)0.063 (-0.636-0.179)0.329 (-0.260-0.539)-1.049 (-1.733-(0.887))
B1-B30.056 (-0.388-0.417)-0.016 (-0.341-0.439)0.391 (-0.093-0.783)-1.140 (-1.452-0.700)
C1.049 (0.775-1.464)1.217 (0.993-1.342)1.627 (1.236-1.718)2.479 (2.266-2.599)
Kruskal-Wallis P<0.001<0.001<0.001<0.001
B1-B3 vs C adjusted P<0.001<0.001<0.001<0.001
B1-B3 vs A-AB adjusted P1.0000.9650.8000.965
C vs A-AB adjusted P<0.001<0.001<0.001<0.001
Masaoka's stage
I-II-0.174 (-0.587-0.188)-0.152 (-0.545-0.132)0.194 (-0.326-0.466)-1.271 (-1.646-(-0.916))
III-IV0.705 (0.314-1.163)0.752 (0.407-1.217)1.096 (0.815-1.599)1.962 (0.707-2.480)
Mann-Whitney P<0.001<0.001<0.001<0.001
Low0.832 (0.552-1.256)1.042 (0.618-1.230)1.316 (0.936-1.633)2.313 (1.909-2.492)
High-0.116 (-0.561-0.189)-0.112 (-0.502-0.179)0.238 (-0.260-0.565)-1.233 (-1.605-(-0.887))
Mann-Whitney P<0.001<0.001<0.001<0.001
Low-0.221 (-0.617-0.155)-0.183 (-0.545-0.063)0.139 (-0.326-0.415)-1.302 (-1.646-(-1.026))
High0.716 (0.475-1.163)0.752 (0.407-1.217)1.096 (0.815-1.599)1.962 (0.693-2.480)
Mann-Whitney P<0.001<0.001<0.001<0.001
Low-0.257 (-0.646-0.111)-0.277 (-0.564-0.011)0.054 (-0.362-0.336)-1.391 (-1.664-(-1.116))
High0.608 (0.314-1.049)0.618 (0.280-1.141)0.936 (0.716-1.485)1.882 (0.643-2.407)
Mann-Whitney P<0.001<0.001<0.001<0.001

Validation of biomarker-related recurrence risk scores

As shown in Supplementary Table 8, all the six biomarkers (miR-130b-5p, miR-1307-3p, miR-425-5p, CD8, CD68, and CCL18) were strongly associated with RFS. However, forward stepwise Cox regression revealed that only miR-425-5p, CD8, and CCL18 were independent prognostic factors for RFS in our cohort. BRRS was significantly greater in WHO type C than in types A-AB or B1-B3 (Table 4) and was positively associated with the risks of recurrence and death as revealed by the univariate Cox regression (HR=2.718, 95% CI: 2.068-3.573, P<0.001; HR=2.779, 95% CI: 1.700-4.544, P<0.001). Multivariate Cox regression showed that BRRS remained significant for RFS but not for OS after being adjusted for gender, age, WHO classification, Masaoka stage, adjuvant chemotherapy, and adjuvant radiotherapy (Table 5 and Figure 4). Time-dependent ROC analysis revealed that BRRS predicted 5-year RFS with better accuracy than miR-425-5p or CCL18 alone (Table 6 and Figure 5). Moreover, the AUC of BRRS for the 5-year RFS was significantly higher than that of the Masaoka stage (AUC: 92.0 vs. 85.2, P=0.007). However, the difference in AUC for OS between CRS and Masaoka stage did not reach significance (Table 6).

Table 5. Multivariate Cox regression for RFS and OS in an independent thymoma cohort.

HR (95% C.I.)PHR (95% C.I.)P
Age1.024 (1.002-1.047).0311.037 (0.997-1.079).069
Gender (Male vs Female)1.387 (0.672-2.862).3770.517 (0.160-1.666).269
Myasthenia gravis (Yes vs No)1.592 (0.788-3.216).1951.353 (0.428-4.274).607
WHO classification (C vs A-B3)0.807 (0.371-1.754).5881.810 (0.437-7.501).413
Masaoka's stage (III-IV vs I-II)14.115 (2.740-72.723).0021.739 (0.087-34.878).718
Adjuvant radiotherapy (Yes vs No)3.450 (1.184-10.058).0230.954 (0.073-12.457).971
Adjuvant chemotherapy (Yes vs No)0.850 (0.235-3.070).8041.34×10-5 (-).979
BRRS2.212 (1.403-3.486).0012.239 (0.886-5.662).089
The BRRS parameters in the table combine multiple biomarkers (including miR-130b-5p, miR-1307-3p, miR-425-5p, CD8, CD68 and CCL18).
Kaplan-Meier curves illustrating the prognostic effects of Masaoka stage (A), WHO classification (B), miR-425-5p (C), CD8 (D), CCL18 (E), and BRRS (F) categorized by corresponding median values. (n = 120 THYM samples).

Figure 4. Kaplan-Meier curves illustrating the prognostic effects of Masaoka stage (A), WHO classification (B), miR-425-5p (C), CD8 (D), CCL18 (E), and BRRS (F) categorized by corresponding median values. (n = 120 THYM samples).

Table 6. Time-dependent ROC analysis showing various accuracies for predicting 3- and 5-year RFS and OS.

RFSTime-dependent AUC (3 year)Time-dependent AUC (5 year)Compared with BRRS at 3 year or 5 yearCompared with CRS at 3 year or 5 year
miR-425-5p0.783 (0.613-0.953)0.875 (0.791-0.960)0.556/0.397
CD80.765 (0.622-0.901)0.832 (0.745-0.920)0.479/0.044
CCL180.737 (0.621-0.852)0.792 (0.709-0.875)0.064/1.183×10-4
BRRS0.815 (0.674-0.956)0.893 (0.816-0.970)0.631/0.275
Masaoka's stage0.776 (0.662-0.890)0.852 (0.780-0.924)0.261/0.007
CRS0.845 (0.748-0.942)0.920 (0.861-0.980)
miR-425-5p0.806 (0.724-0.887)0.790 (0.592-0.987)0.334/0.415
CD80.693 (0.421-0.965)0.710 (0.509-0.911)0.507/0.298
CCL180.781 (0.731-0.831)0.697 (0.521-0.874)0.951/0.067
BRRS0.771 (0.678-0.863)0.776 (0.580-0.972)0.686/0.703
Masaoka's stage0.813 (0.764-0.861)0.726 (0.550-0.902)0.465/0.349
CRS0.774 (0.684-0.864)0.773 (0.573-0.973)
Time-dependent ROC illustrating the accuracy of Masaoka stage, BRRS, and CRS for predicting the RFS at 3 years (A) and 5 years (B). (n = 120 THYM samples).

Figure 5. Time-dependent ROC illustrating the accuracy of Masaoka stage, BRRS, and CRS for predicting the RFS at 3 years (A) and 5 years (B). (n = 120 THYM samples).


TETs arise from thymic epithelial cells with a low morbidity rate. However, few studies have been done on the treatment and prognosis of those patients presenting with TETs. Here, we developed BRRS for TETs using microRNA-mRNA pairs and tumor-infiltrating immune cells from the TCGA dataset. Our analysis of 97 tumor samples from the TCGA thymoma datasets indicated that macrophage infiltration abundance was the most pronounced prognostic factor for PFI. A total of 41 microRNA-mRNA pairs were further analyzed to obtain 4 microRNAs responsible for biomarkers of PFI. Then, we analyzed the clinical features and found that the prognostic indicators consisting of microRNAs were superior to the current clinical staging.

IMRS3 and CD8+ T cells were found related to the TREM2hi score. In addition, we validated miR-435-5p, CD8, and CCL18 as independent predictors of RFS. In the 99 cases, BRRS was significantly related to the recurrence risk and death risk. Notably, after adjusting for gender, age, myasthenia gravis, Masaoka staging, WHO classification, radiotherapy, and chemotherapy, BRRS remained an independent prognostic factor for RFS. The tumor microenvironment has important interactions with cancer cells to determine outcomes, and especially immune-related microRNAs are important biomarkers for prognosis and treatment response [11]. The immune-related microRNA signature was strongly correlated with the immune microenvironment and clinical outcomes of hepatocellular carcinoma [12]. The expression of immune-related microRNAs was associated with the role of the tumor immune microenvironment in immunosuppressive therapy [13]. In this study, BRRS, independent of Masaoka staging, had a significant regression contribution to RFS, suggesting that BRRS may serve as a promising signature in addition to the current clinical prognostic indicators.

Additionally, we validated that IMRS3 is related to TREM2hi and acts as a prognostic factor superior to the current clinical staging. Besides, IMRS of macrophages was significantly correlated with the thymoma PFI collected from the TCGA dataset. Immune-related genes and microRNAs have been widely studied, and they hold significate research value for cancer prognosis. Immune-related genes can effectively predict the clinical outcome of hepatocellular carcinoma treatment and assess the clinical status of the tumor [14]. Guo et al. developed an immune-based prognostic model for renal cell carcinoma, validating that the immune-related microRNA signature is a feasible and reliable prognostic model [15]. Immune-related microRNAs are non-invasive biomarkers of clinical outcomes and modulate the immunosuppressive tumor microenvironment in metastatic melanoma [11]. Immune-related microRNAs regulate immune cells directly [15] or may modulate the immune microenvironment through exosomes. Unlike the study in which prognostic markers were constructed with potential microRNAs, the present study applied a more indirect statistical-based approach to screen microRNAs that may modulate macrophage function. We successfully validated that miR-130b-5p, miR-1307-3p, and miR-425-5p in tumors were significantly related to the WHO classification, and their increased expression levels may be associated with the development of immunosuppressive microenvironment in thymoma. All IMRS, WHO classification, and Masaoka staging were prognostic factors for RFS and OS. Therefore, we believe that immune-related microRNAs and IMRS may be vital contributors in predicting the prognosis of thymoma. In this sense, the immune-related microRNAs we identified differ from the “immune-related microRNAs” defined by microRNA species that directly affect gene regulation in immune cells.

In the 120 thymomas, IMRS3 and TREM2hi scores were found significantly positively correlated, while TREM2hi was negatively correlated with CD8+ T lymphocytes. Thymoma led to polarization of M2 macrophages due to high expression of IMRS3, which rejected CD8 T and other effector lymphocytes and promoted the malignant progression of thymoma. The high TREM2hi score expression group was mainly characterized by immunoregulatory pathways such as the IL-10 signaling pathway and IFN-γ pathway. The expression of hsa-microRNA-1307-3p and hsa-miR-425-5p, which were significantly correlated with PFI in IMRS3, had a strong positive correlation with CCL18. Hence, we proposed that immune-related microRNAs may be a biomarker for the immunosuppressive tumor microenvironment through modulating IL-10 or IFN-γ-mediated signaling pathway. Immune-related microRNAs regulate immunoregulatory molecules, which mediate immune escape from cancer and result in a poor prognosis. Immunoregulatory microRNAs, as a novel mechanism, can reverse immune escape from tumors [16]. Proinflammatory cytokines (interleukin and tumor necrosis factor-mediated signaling pathways) are regulated by the expression of immune-related microRNAs and associated with colorectal tumor progression, which can be considered as biomarkers of colorectal cancer (immune microenvironment of colorectal tumors: involvement of immune genes and microRNAs belonging to the TH17 pathway). Moreover, IMRS3 expression was significantly associated with CD8 expression in the IHC data of our clinical samples. Therefore, we believe that immune-related microRNAs can act as an accurate prognostic marker for thymoma in the future.

CD8 is a marker of killer T lymphocytes and CD68 is a marker of tumour-associated macrophages (TAM). CCL18 is primarily produced and secreted by the innate immune system (including TAM) and affects primarily the adaptive immune system (lymphocytes). These three interact to ultimately recruit Th2 and regulatory T cells, which mediate immunosuppressive effects; and secrete a variety of cytokines, growth factors (e.g. EGF, VEGF, PDGF, FGF and TGFβ), which interact with and influence multiple cell types in the tumor microenvironment and, in turn, mediate complex effects. Thus, the triad of CD8, CD68 and CCL18 clearly influences tumor prognosis, which corroborates our findings in the cox regression analysis.

Our study still has some limitations. First, we used the MCPcounter method to analyze immune cell abundance, which simply utilized the geometric means of the gene tag expression values of 10 immune cells to indirectly measure the abundance of infiltration, which may not represent the physical meaning of the actual infiltration ratio. Second, although the immune-related RNAs we proposed show great predictive value for the prognosis of thymoma, more in-depth studies are highly warranted to validate this finding. In addition, the immune-related microRNAs we identified in the present study need to use scRNAseq data to determine in a more accurate way whether expression of thymoma epithelial cells and microenvironmental macrophages via exosomes and other methods lead to an immunosuppressive environment. The specific mechanism needs to be further investigated. Despite these shortcomings, our research still has important scientific implications. For the first time, we propose that immune-related microRNAs play a critical role in thymoma prognosis and have great potential for clinical use.

Immune-related microRNAs may serve as essential biomarkers for predicting the prognosis of thymoma, and they may affect the prognosis by modulating the immune microenvironment. Our efforts provide a promising approach for accurate judgment of thymoma progression, leading to a new research interest that immune-related microRNAs that mediate the immune microenvironment may contribute to poor prognosis and even pathogenesis of thymoma.

Materials and Methods

Dataset retrieval and preprocessing

Gene expression values, as well as phenotype files and survival data for 122 thymomas, were downloaded from UCSC (, and the expression matrices were in a format of log2(x+1)-transformed RSEM normalized counts. Genes with 0 count more than 20% of the samples were filtered, and 16,552 genes were included. Additionally, we also downloaded microRNAs expression matrices for 126 thymomas from UCSC in log2(RPM+1) format, and 411 microRNAs were available for analysis after filtering the microRNAs with NA-containing expression values. Accession ID numbers of microRNAs were transformed with mature microRNA names using the miRBaseConverter package [17], and a total of 408 annotated microRNAs were finally included.

Evaluation of the abundance of tumor-infiltrating immune cells

The abundance of 8, 22, and 10 immune cells was analyzed using the EPIC, CIBERSORTx, and MCPcounter methods, respectively [18, 19]. For EPIC, count values were used as imported data to obtain the abundance of 8 cell populations (including CD4 T cells, CD8 T cells, NK cells, and macrophages) via deconvolution. Similar to EPIC, CIBERSORTx also utilized deconvolution to obtain 22 immune cell abundances. In the MCPcounter process, the infiltration abundance of 10 immune cells was indirectly measured with the geometric mean of their tag expression values.

Identification of microRNA-mRNA pairs influencing macrophage function on relapse-free survival

In order to screen out the microRNAs and gene mRNAs that may indirectly affect RFS by affecting tumor immune cell function, a similar method for constructing T cell depletion signatures proposed by Jiang [20] was used. Subsequently, a prognostic microRNA-mRNA regulatory network for the identification method was constructed using the microRNA-mRNA target genes.

By univariate Cox regression, “macrophage” abundance obtained by EPIC was found to be a significant prognostic factor for PFI (β = 90.08, P = 0.0113). Given this, regression analysis of the 408 microRNAs was performed using macrophage abundance, microRNA expression values, and Cox regression of interaction terms between the two, and mRNA expression values of the 16,552 genes underwent Cox regression in a similar manner. MicroRNA species and mRNA genes with macrophage abundance, expression values, and interaction term regression coefficients Wald P < 0.05 were extracted as candidate genes.

The correlation of the expression values of these candidate genes was further analyzed with the R package anamiR function “negativecor” using the Pearson correlation coefficient method. MicroRNA-mRNA target gene pairs with r < −0.5 were included in the “databasesupport” function [21] to obtain microRNA-mRNA target gene pairs predicted by the miRNA database, and the microRNAs whose interaction with mRNAs was supported by at least three datasets were finally selected as candidates for weighted linear modeling. Penalized Cox regression was performed using LASSO to obtain a PFI prognostic signature containing four microRNA species, which was named the immune-related microRNAs score (IMRS). In order to make validation more feasible in a clinical setting, a reduced model incorporating the first three microRNAs (IMRS3) was calculated in comparison to the entire model (IMRS4) in terms of PFI prognostic performance. Fivefold cross-validation was conducted to select the lamda using the function cv.glmnet [22]. IMRS was calculated by summing the expression levels of microRNAs weighted with the corresponding non-zero coefficients.

Characterization of the immune microenvironment highly infiltrated by M2 macrophages

To further reveal the correlation of IMRS with M2 macrophages and T cells, a gene signature peculiar to TREM2hi macrophages was utilized [23], which consists of 11 genes (TREM2, SPP1, RNASE1, MT1G, SEPP1, FOLR2, NUPR1, KLHDC8B, CCL18, MMP12, and APOC2). The mean expression value of the 11 genes was taken as the TREM2hi score.

The correlation between the 11 genes, TREM2hi score, infiltrating immune cell abundance, and IMRS was analyzed using the Pearson method. In order to further explore the effect of TREM2hi score on tumor immune microenvironment, 120 thymoma samples were divided by the median TREM2hi score into high and low expression groups, and differentially expressed genes were obtained using “limma” package “voom” and “lmFit”. The enrichment analysis of the Reactome gene set was done using the clusterProfiler package “gesPathway” function with the t as the test variable [24], and a corrected P < 0.05 was taken as the significantly enriched pathway.


To explore the correlation between the expression of the three candidate microRNAs, immune cells, and RFS, 99 patients with histologically confirmed thymoma or thymic carcinoma who were admitted into the Department of Thoracic Surgery at Daping Hospital of the Third Military Medical University between January 2013 and April 2017 were enrolled in this study. Those with malignancies were excluded. The tumors were classified according to the 2015 WHO criteria and then staged as per the modified Masaoka staging system. Approval for this study was gained from the Medical Ethics Committee at the hospital (approval no. 2021-100), and written informed consent was obtained from all patients prior to enrolment. Patient information, including age, gender, histological findings, presence of myasthenia gravis, Masaoka stage, WHO histological type, and follow-up information, was collected from clinical records and by questionnaires.

Quantitative RT-PCR (qRT-PCR)

We extracted total RNA from the 99 formalin-fixed, paraffin-embedded specimens with the Tiangen RNAprep Pure FFPE Kit (TIANGEN, Beijing, China, Category Number: DP439), detected the RNA concentration and purity with Nano Drop® 2000 (Thermo Fisher Scientific, CA, USA), reversely transcribed cDNA with the HiScript 1st Strand cDNA Synthesis Kit, and then conducted quantitative PCR reactions using AceQ qPCR SYBR Green Master Mix (without ROX) (Vazyme Biotech Co., Ltd, Nanjing, China) and Bio-Rad CFX96 system (Bio-Rad Laboratories, Inc., China). Each sample was assayed in triplicate, and U6 snRNA was used for normalization. The expression level of individual microRNA was determined by the 2-ΔΔCT method (ΔCT = CTmicroRNA – CTU6RNA), and one patient specimen was randomly selected as a reference sample. The primers used in this study are shown in Supplementary Table 1.


For immunohistochemistry analysis, the paraffin-embedded tissue blocks were cut into sections at 3 μm thickness, and human CCL18 polyclonal antibodies (catalog no. ab233099, Cambridge, MA, USA) as well as CD68 and CD8 primary antibodies (catalog nos. ZM-0060 and ZA-0508, Zhong Shan Golden Bridge Biological Technology, Beijing, China) were used. Subsequently, immunostaining was done with the SPlink Detection kit (Biotin-Streptavidin HRP Detection System; catalog no. SP-9001; Zhong Shan Golden Bridge Biological Technology, Beijing, China) according to the manufacturer’s protocol. Slides were examined by two independent pathologists (Department of Pathology, Daping Hospital, Army Medical University) who had no prior knowledge of the clinical and pathological parameters. The staining intensity was classified as follows: grade 0: no staining; grade 1: definite but weak staining; grade 2: moderate staining; grade 3: strong staining. Positive cells was counted at ×400 magnification (range, 0-100%) in five random microscopic fields and categorized as follows: grade 0: no staining; grade 1: ≤10%; grade 2: 11%-50%; grade 3: 51%-75%; grade 4: >75%. The final scores were calculated by multiplying the proportion score by the intensity score, which were classified into 4 grades: 0, negative (-); 1-4, (+); 5-8 (++); 9-12 (+++). Samples that were ++ or +++ were defined as highly expressed and the rest as lowly expressed.

Statistical analysis

All microRNA expression levels were transformed logarithmically and expressed as median values and interquartile ranges. The differences in microRNA expression among the WHO classifications, between Masaoka stage I-II and Masaoka stage III-IV subgroups, or between high and low expression subgroups of CD8, CD68, and CCL18 were evaluated by Kruskal-Wallis test or Mann-Whitney U test. The χ2 test or Fisher’s exact probability was applied to analyze the correlation between CD8, CD68, CCL18, and clinicopathological factors. The RFS of different subgroups was compared using the Kaplan-Meier method and the log-rank test, and prognostic factors for RFS and OS were identified in univariate and multivariate Cox proportional hazard regression models. In particular, a likelihood ratio test was used to determine whether a covariate was significantly entered into the regression model. Two weighted linear models were developed, one of which was based on the biomarkers in question (miR-130b-5p, miR-1307-3p, miR-425-5p, CD8, CD68, and CCL18), i.e., BRRS, and the other was based on a combination of BRRS and clinical features, i.e., combined risk scores (CRS). The predictive efficiency of Masaoka staging, BRRS, CRS, and other individual biomarkers for 3- and 5-year RFS and OS was determined with time-dependent ROC curve analysis by using the function “timeROC” and the function “compare” implemented in the R package. All other statistical analyses were performed using the SPSS 17.0 software (IBM SPSS, Chicago, IL, USA). All tests were bilateral, and P < 0.05 was considered statistically significant.

Ethical disclosure

The study was authorized by the ethics Committee of Army Medical Center of PLA Approval of Medical Research Involving People Ethical Ratification NO:2021 (100).


TETs: Thymic epithelial tumors; WHO: World Health Organization; qRT-PCR: Quantitative RT-PCR; CRS: combined risk scores.

Author Contributions

All authors had full access to the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. Study concept and design: H. X, and B.W. Acquisition of data: R. Y, H. C, X. Y, B.W and Z. Y. Performed experiments: B. W, H. X, C. C, JX. C, and ZM. Z. Analysis and interpretation of the data: H. X, B.W Drafting of the manuscript: B. W, H. X, C. C, JX. C, and ZM. Z. Critical revision of the manuscript for important intellectual content: B. W, H. X, C. C, JX. C, and ZM. Z. Statistical analysis: B. W, H. X, C. C, JX. C, and ZM. Z. Obtained funding: B.W. Administrative, technical and material support: B.W and B.W. Study supervision: B.W.


We would like to acknowledge the reviewers for their helpful comments on this paper.

Conflicts of Interest

The authors declared that they have no conflicts of interest to this work.

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.


This work was supported by National Natural Science Foundation of China (No. 81672536,82173170) Science and Technology Research Program of Chongqing Municipal Education Commission (No. KJZD-K202000405), Innovation research group in Colleges and Universities Program of Chongqing Municipal Education Commission (No. CXQT20012) and the project of the top-notch talent cultivation program for the graduate students of Chongqing Medical University No.BJRC201918.


  • 1. Engels EA. Epidemiology of thymoma and associated malignancies. J Thorac Oncol. 2010 (Suppl 4); 5:S260–5. [PubMed]
  • 2. Kelly RJ, Petrini I, Rajan A, Wang Y, Giaccone G. Thymic malignancies: from clinical management to targeted therapies. J Clin Oncol. 2011; 29:4820–7. [PubMed]
  • 3. Travis WD, Brambilla E, Burke AP, Marx A, Nicholson AG. Introduction to The 2015 World Health Organization Classification of Tumors of the Lung, Pleura, Thymus, and Heart. J Thorac Oncol. 2015; 10:1240–2. [PubMed]
  • 4. Kondo K, Yoshizawa K, Tsuyuguchi M, Kimura S, Sumitomo M, Morita J, Miyoshi T, Sakiyama S, Mukai K, Monden Y. WHO histologic classification is a prognostic indicator in thymoma. Ann Thorac Surg. 2004; 77:1183–8. [PubMed]
  • 5. Feng Y, Lei Y, Wu X, Huang Y, Rao H, Zhang Y, Wang F. GTF2I mutation frequently occurs in more indolent thymic epithelial tumors and predicts better prognosis. Lung Cancer. 2017; 110:48–52. [PubMed]
  • 6. Qadir MI, Faheem A. miRNA: A Diagnostic and Therapeutic Tool for Pancreatic Cancer. Crit Rev Eukaryot Gene Expr. 2017; 27:197–204. [PubMed]
  • 7. Chen C, Liu JM, Luo YP. MicroRNAs in tumor immunity: functional regulation in tumor-associated macrophages. J Zhejiang Univ Sci B. 2020; 21:12–28. [PubMed]
  • 8. O’Donnell JS, Teng MW, Smyth MJ. Cancer immunoediting and resistance to T cell-based immunotherapy. Nat Rev Clin Oncol. 2019; 16:151–67. [PubMed]
  • 9. Korbecki J, Olbromski M, Dzięgiel P. CCL18 in the Progression of Cancer. Int J Mol Sci. 2020; 21:7955. [PubMed]
  • 10. Manjili MH. Tumor Dormancy and Relapse: From a Natural Byproduct of Evolution to a Disease State. Cancer Res. 2017; 77:2564–9. [PubMed]
  • 11. Jorge NA, Cruz JG, Pretti MA, Bonamino MH, Possik PA, Boroni M. Poor clinical outcome in metastatic melanoma is associated with a microRNA-modulated immunosuppressive tumor microenvironment. J Transl Med. 2020; 18:56. [PubMed]
  • 12. Li Y, He X, Zhang X, Xu Y, Wu Y, Xu X. Immune-related microRNA signature for predicting prognosis and the immune microenvironment in hepatocellular carcinoma. Life Sci. 2021; 265:118799. [PubMed]
  • 13. Tokumaru Y, Asaoka M, Oshi M, Katsuta E, Yan L, Narayanan S, Sugito N, Matsuhashi N, Futamura M, Akao Y, Yoshida K, Takabe K. High Expression of microRNA-143 is Associated with Favorable Tumor Immune Microenvironment and Better Survival in Estrogen Receptor Positive Breast Cancer. Int J Mol Sci. 2020; 21:3213. [PubMed]
  • 14. He Y, Dang Q, Li J, Zhang Q, Yu X, Xue M, Guo W. Prediction of hepatocellular carcinoma prognosis based on expression of an immune-related gene set. Aging (Albany NY). 2020; 12:965–77. [PubMed]
  • 15. Guo Y, Li X, Zheng J, Fang J, Pan G, Chen Z. Identification of a novel immune-related microRNA prognostic model in clear cell renal cell carcinoma. Transl Androl Urol. 2021; 10:888–99. [PubMed]
  • 16. Yi M, Xu L, Jiao Y, Luo S, Li A, Wu K. The role of cancer-derived microRNAs in cancer immune escape. J Hematol Oncol. 2020; 13:25. [PubMed]
  • 17. Xu T, Su N, Liu L, Zhang J, Wang H, Zhang W, Gui J, Yu K, Li J, Le TD. miRBaseConverter: an R/Bioconductor package for converting and retrieving miRNA name, accession, sequence and family information in different versions of miRBase. BMC Bioinformatics. 2018 (Suppl 19); 19:514. [PubMed]
  • 18. Racle J, de Jonge K, Baumgaertner P, Speiser DE, Gfeller D. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife. 2017; 6:e26476. [PubMed]
  • 19. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, Selves J, Laurent-Puig P, Sautès-Fridman C, Fridman WH, de Reyniès A. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016; 17:218. [PubMed]
  • 20. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, Liu J, Freeman GJ, Brown MA, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018; 24:1550–8. [PubMed]
  • 21. Wang TT, Lee CY, Lai LC, Tsai MH, Lu TP, Chuang EY. anamiR: integrated analysis of MicroRNA and gene expression profiling. BMC Bioinformatics. 2019; 20:239. [PubMed]
  • 22. Sveen A, Ågesen TH, Nesbakken A, Meling GI, Rognum TO, Liestøl K, Skotheim RI, Lothe RA. ColoGuidePro: a prognostic 7-gene expression signature for stage III colorectal cancer patients. Clin Cancer Res. 2012; 18:6001–10. [PubMed]
  • 23. Xiong D, Wang Y, You M. A gene expression signature of TREM2hi macrophages and γδ T cells predicts immunotherapy response. Nat Commun. 2020; 11:5084. [PubMed]
  • 24. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. [PubMed]