A novel immune-related microRNA signature for prognosis of thymoma

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 TREM2^hi 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.

Characteristics of immune microenvironment in thymoma

As shown in Figure 1, there was a clear trend that IMRS3 and the TREM2^hi 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 TREM2^hi signature (Figure 2). In addition, CCL18 (Chemokine (C-C motif) ligand 18), one of the components of the TREM2^hi signature, was also highly positively correlated with miR-1307-3p, miR-425-5p, IMRS3, and TREM2^hi signature itself (Figure 2). Gene set enrichment analysis showed that tumors with high TREM2^hi 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 TREM2^hi macrophages, may be essential for thymoma progression.

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.

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 TREM2^hi 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.

Figure 3. 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).

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

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

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

Dataset retrieval and preprocessing

Gene expression values, as well as phenotype files and survival data for 122 thymomas, were downloaded from UCSC (https://xenabrowser.net/datapages/), 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 TREM2^hi 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 TREM2^hi score.

The correlation between the 11 genes, TREM2^hi score, infiltrating immune cell abundance, and IMRS was analyzed using the Pearson method. In order to further explore the effect of TREM2^hi score on tumor immune microenvironment, 120 thymoma samples were divided by the median TREM2^hi 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.

Patients

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 = CT_microRNA – CT_U6RNA), and one patient specimen was randomly selected as a reference sample. The primers used in this study are shown in Supplementary Table 1.

Immunohistochemistry

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 χ² 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).