Associated with Immune Function, miR-150-5p is a Favorable Biomarker for Head and Neck Cancer Patients

Gene


Abstract
Background Accumulating evidence has shown that dysregulated expression of microRNAs plays a key role in tumorigenesis. To explore the mechanisms of this we conducted this study.

Methods
Five Gene Expression Omnibus datasets (GEO) datasets , GSE32960, GSE36682, GSE43039, GSE70970 and GSE118613 and head and neck squamous cell carcinoma data of The Cancer Genome Atlas (TCGA) were analysis in this study.

Results
By analyzing the microRNA expression pro le of nasopharyngeal carcinoma (NPC) in the ve GEO datasets, we identi ed miR-150-5p as potential biomarker for patient survival. To explore the mechanisms of this, We examined the head and neck squamous cell carcinoma data of TCGA and found that miR-150-5p was correlated with high enrichment of tumor-in ltrating B cells, low enrichment of cancer-associated broblasts and down-regulated oncogenic pathways. miR-150-5p may also improve the immune response in the tumor microenvironment. These ndings may explain how miR-150-5p improves outcome of head and neck squamous cell carcinoma patients including NPC. Additionally, the exosomal long noncoding RNA AC073130.1 was identi ed as a potential regulator of miR-150-5p. As miR-150-5p can also be released via exosomes, this study provides insight into the cross-talk of tumor cells and B cells in the tumor microenvironment via exosomal AC073130.1 and miR-150-5p.

Conclusion
MiR-150-5p improves outcome of head and neck squamous cell carcinoma patients by improving the immune response. There might be a cross-talk of tumor cells and B cells in the tumor microenvironment via exosomal AC073130.1 and miR-150-5p.

Background
Nasopharyngeal carcinoma (NPC), prevailing in the Southeast Asia, is an epithelial carcinoma arising from nasopharyngeal mucosa [1]. Although signi cant advances have been achieved in the treatment strategies for NPC, NPC is a major cause of death in epidemiological areas [1].
MicroRNAs (miRNAs) are a type of endogenous non-coding RNA approximately 21-23 nucleotides in length. miRNA can inhibit the translation of target mRNAs by forming a miRNA-induced silencing complex [2]. Studies have shown that miRNA can function as tumor suppressors or oncogenes in many types of cancer including NPC by affecting cell cycle, proliferation, apoptosis, migration, invasion and radiotherapy sensitivity [3,4] .
To gain more insights into the mechanism underlying NPC, we investigated the miRNA expression pro le in NPC by analyzing ve Gene Expression Omnibus (GEO) datasets. miR-150-5p was found as tumor suppressor and was down-regulated in NPC tissues compared with non-NPC tissues. However the function of miR-150-5p is controversial. For example, miR-150-5p was found to exhibit both tumor suppressor functions [5] as well as oncogenic functions in non-small cell lung cancer [6]. miR-150-5p was reported as an exosomal non-coding RNA [7] and can be released by components in tumor microenvironment (TME) [8]. Therefor we hypothesis that miR-150-5p can not only affect the biological behavior of tumor cells, but also affect TME around tumor cells. Then we analyzed samples in head and neck squamous cell carcinoma (HNSC) data of The Cancer Genome Atlas (TCGA) database to explore the potential anti-cancer mechanism of miR-150-5p.

Datasets
Five datasets, GSE32960, GSE36682, GSE43039, GSE70970 and GSE118613, including miRNA-seq data of NPC in GEO, were analyzed. GSE118613 data were from blood samples. Survival status, survival time, metastasis status, and time to metastasis in data from GSE70970 were used for analyses. miRNA-seq data and mRNA-seq data of HNSC patients in TCGA who had miRNA-seq data available at the same time were collected. The clinical characteristics of the HNSC patients were also collected.

Bioinformatics analysis
Analyses were performed using R software 3.6.3. DESeq2 was used for differentially expressed gene (DEG) analysis. miRNAs with LogFC > 1 and adjusted p < 0.1 were selected as miRNAs with signi cant expression [9]. For mRNAs, LogFC > 1 and adjust p < 0.05 were considered to indicate mRNAs with signi cant expression. Volcano plot and heatmap of differentially expressed genes were obtained by ggplot2.
Using the median expression of miR-150-5p, HNSC patients were divided into the high miR-150-5p expression group and low miR-150-5p expression group. Survival package and survminer package were used for Kaplan-Meier survival plots and log-rank test was used to compare the survival curves. P < 0.05 was considered to indicate statistical signi cance. Cox regression analysis was performed using the survival package.
Gene Ontology (GO) enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Set Enrichment Analysis (GSEA) were performed using clusterPro ler package.
Correlation analysis was tested by Spearman correlation.

Tumor-in ltrating cells
The tumor-in ltrating cells in HNSC were calculated by Timer [10], quanTIseq, MCP-Counter, EPIC, xCell [11] and CIBERSORT [12] respectively. The enrichment of different cells was compared by Wilcox test. The in ltrating cells that were signi cantly up-regulated or down-regulated were selected. GSVA package [13] was used for the correlation of AC073130.1 and B cells.

Network of ceRNAs
The lncRNAs that regulate miR-150-5p were identi ed by rst determining the predicted miR-150-5pregulating lncRNAs using LncBase [14] and the highly expressed lncRNAs in patients with low miR-150-5p expression and then identifying the overlapping lncRNAs. The potential mRNAs regulated by miR-150-5p were identi ed by rst determining the predicted mRNA targets of miR-150-5p by TargetScan [15] and the highly expressed mRNA in patients with low miR-150-5p expression and then identifying the overlapping mRNAs.

miR-150-5p showed decreased expression in NPC patients and was associated with poor prognosis
To examine the miRNA expression pro le in NPC, we performed analysis of ve NPC datasets in GEO, GSE32960, GSE36682, GSE43039, GSE70970 and GSE118613.Two microRNAs, miR-150-5p and miR-342-3p, were consistently down-regulated in the NPC group compared with the non-NPC group in all ve GEO datasets ( Fig. 1A-B, Fig. S1A). Decreased expression of miR-150-5p and miR-342-3p was associated with poor overall survival (OS) in NPC patients (Fig. 1C, Fig. S1B). Decreased miR-150-5p expression was also associated with poor metastasis-free survival in NPC patients (Fig. 1D). miR-342-3p had no effect on the metastasis-free survival in NPC patients (Fig. S1C).
Gene expression pro le associated with miR-150-5p We next explored the potential mechanism by which miR-150-5p affects NPC patient survival using TCGA HNSC database. miR-150-5p was decreased in HNSC tissues compared with paired non-tumor tissues ( Fig. 2A). Decreased expression of miR-150-5p associated with poor OS and progression-free survival (PFS) in HNSC patients ( Fig. 2B-C). miR-150-5p was also an independent prognostic factor of OS and PFS in HNSC (Fig. 2D).
We next compared the gene expression pro les of HNSC with high miR-150-5p expression and those with low miR-150-5p expression (Fig. 2E-G). GO enrichment analysis of the differentially expressed genes were mainly associated with immune-related terms, such as lymphocyte differentiation, T cell activation, B cell activation, lymphocyte-mediated immunity, humoral immune response, antigen-receptor mediated signal pathway, B cell-mediated immunity, B cell proliferation, and B cell receptor signal pathway ( Fig. 2H and Supplementary table 1). KEGG pathway analysis revealed that the differentially expressed genes were also mainly involved in immune-related pathways, such as natural killer cell-mediated cytotoxicity, primary immunode ciency, T cell receptor pathway, antigen processing and presentation, Th17 cell differentiation, B cell receptor signal pathway, graft versus host disease, intestinal immune network for IgA production, and allograft rejection ( Fig. 2I and Supplementary table 1).
GSEA analysis showed increased activated pathways in B cell activation, B cell-mediated immunity, leukocyte-mediated immunity, activation of immunity, immune response to tumor cell, allograft rejection, graft versus host disease, B cell receptor signal, antigen processing and presentation, T cell receptor signal, complement cascade and adaptive immune system in HNSC with high expression of miR-150. While the activation of pathways in MET signal, cell junction organization, non-integrin membrane-ECM interactions, hypoxia, glycolysis, angiogenesis and epithelial mesenchymal transition were decreased in HNSC with high expression of miR-150. (Fig. 2J-N (Fig. 3A). Exosomes function as messengers of cell-cell communication and are released into paracellular and body uids, such as plasma, salvia, urine, and breast milk [20,21]. The cross-talk between tumor cells and stromal cells via exosomes in the TME regulates TME remodeling, invasion, metastasis, angiogenesis, immune suppression and drug resistance [22,23]. Therefore, we investigate the correlation of miR-150-5p and TME.
The proportions of TME-in ltrating cells in HNSC were calculated using Timer, quanTIseq, MCP-Counter, CIBERSORT and xCell were used to calculate the abundance of the subtypes of B lymphocytes. Naïve B cells, memory B cells and plasma B cells were all positively correlated with miR-150-5p expression (Supplementary table 3). Higher level of plasma B cells was also correlated with longer survival of HNSC patients ( Fig. 3N and Supplementary table 4). These results suggested that the impaired survival of HNSC patients may be related to the decrease of B lymphocyte in the TME. miR-150-5p was regulated by AC073103.1 As miR-150-5p was found down-regulated in NPC and TCGA HNSC compared with non-tumor tissues, we next investigate the potential cause of decreased miR-150-5p. To explore potential ceRNAs that may regulate miR-150-5p, we examined lncRNA expressions in HNSC with high miR-150-5p expression and those with low miR-150-5p expression (Fig. 4A-C). Subsequent analysis identi ed 89 potential lncRNAs that regulate miR-150-5p and 6 mRNAs that may be regulated by miR-150-5p (Fig. 4D, E). The ceRNA network associated with miR-150-5p is shown in Fig. 4F.
The exoRBase database is a repository of circular RNA (circRNA), lncRNA and mRNA derived from RNAseq data analyses of human blood exosomes [24]. The lncRNA AC073103.1 was one of the candidate regulators of miR-150-5p identi ed above and was found in blood exosomes in the exoRBase database (Fig. 4G). AC073103.1 was increased in cancer patients compared with normal person (Fig. 4G). In TCGA HNSC, the expression of AC073103.1 was increased in tumor tissues compared with paired non-tumor tissues (Fig. 4H). High expression of AC073103.1 was associated with poor survival in HNSC. Furthermore, high expression of AC073103.1 was associated with low enrichment of B cells (Fig. 4J-L). These data indicated that tumor cells may down-regulated immune reaction pathways and decrease the enrichment of immune-related cells in the TME by increasing the expression of AC073103.1 (Fig. 4M) .

Discussion
Dysregulation of miRNAs plays a key role in cancer progression in most cancers including NPC [25]. In addition to human miRNAs, miRNAs encoded by Epstein-Barr virus (EBV) also affect the proliferation, invasion, and immune evasion of tumor cells [25]. Only GSE43039[26] contained EBV-encoded miRNAs among all the ve GEO datasets, and therefore no EBV-encoded miRNAs were analyzed in this study. miR-150-5p was down-regulated in GEO NPC and TCGA HNSC datasets. miR-150-5p was also associated with OS and DFS of NPC as well as OS and PFS of HNSC. Our ndings suggest that a higher immune response, lower enrichment of CAF and down-regulation of oncogenic pathways may represent mechanisms by which high expression of miR-150-5p led to better patient prognosis.
In this study most of the differentially expressed genes between high and low miR-150-5p expression groups were enriched in immune-related GO terms and KEGG pathways. GSEA analysis also showed that the immune response was higher in patients with higher miR-150-5p expression. As immune suppression is one of the main factors that cause cancer immune evasion and tumor development [27], a suppressed immune response may be the most important potential cause of poor outcomes in patients with low miR-150-5p expression.
Recent studies has revealed that tumor-in ltrating B cells could promote the T cell response by forming a tertiary lymphoid structure[28] and improve immunotherapy response and survival [29]. In this study, many B cell-related pathways were activated in the miR-150-5p high expression group and B cells were correlated with the expression of miR-150-5p in all ve predicting algorithms. B cells were also associated with the survival of HNSC patients. Furthermore, plasma B cells, the effector cell of B cells, also correlated with the expression of miR-150-5p and survival of HNSC patients. This evidence suggested that tumorin ltrating B cells in HNSC were responsible for patient outcomes.
Previous studies have indicated that miR-150 may function to maintain homeostasis of B cell development: as an exosomal miRNA, miR-150-5p is highly expressed in mature B cells, while it blocks maturation of premature B cells [30,31]. Therefore, high level of miR-150-5p in TME could not increase the enrichment of mature B cells, but it could illustrate increased enrichment of mature B cells. However, some studies showed that high expression of miR-150 enhances the antigen presentation and B cell receptor signaling [32,33]. Therefore high expression of miR-150-5p might up-regulate the immune response and improve patient outcome.
In this study, a lower abundance of CAFs was found in the high miR-150-5p expression group. CAFs originate from normal broblasts, and hypoxia-inducible factor (HIF)-1α transcription factor is one of the main factors that promote the conversion of normal broblasts to CAFs [34]. The down-regulated pathway of hypoxia in this study might be the reason underlying the low abundance of CAFs in the high miR-150-5p expression group. Recent studies have shown the tumor-promoting roles of CAFs, such as functions in promoting the proliferation and stemness of cancer, facilitating invasion and migration, controlling vascular and immune system and drug resistance [17,35]. Therefore, a reduced abundance of CAFs might be another reason of better outcome in patients with high miR-150-5p expression.
The GSEA analysis also revealed that several oncogenic pathways were down-regulated in the case of high expression of miR-150. Although miR-150-5p acts as oncogene in some cancers[6, 36-38], miR-150-5p has mostly been show to play an anti-tumor role in cancer [5,[39][40][41][42][43][44][45][46][47]. The different roles of miR-150-5p in different studies and cancer types might depend on the target gene it regulates. The anti-cancer role of miR-150-5p in our study is consistent with others in HNSC [40,42,45]. Some lncRNAs contain binding sequences for speci c target miRNAs and thus can act as 'sponges' or competitive endogenous RNAs to affect the subsequent activation of target mRNAs [48]. Studies have reported that a large number of RNA transcripts contain miRNA-binding sites, indicating that these RNA transcripts may regulate each other by competing for shared target miRNAs [49]. Additionally, dysregulation of non-coding RNAs is often involved with cancer [50]. Therefore, we explored the potential lncRNAs that may sponge miR-150-5p. Among the overlapping down-regulated lncRNAs in the miR-150-5p high expression group and lncRNAs predicted by LncBase, one lncRNA, AC073103.1, was correlated with patient survival and enrichment of B cells in the TME. Previous studies showed that AC073103.1 is released via exosomes and up-regulated in some types of cancers [24]. Therefore, tumor cells may reshape the TME by releasing exosomal AC073103.1. We also identi ed six potential target mRNAs of miR-150-5p and three, GHSR [51], IGF2BP1 [52], and PEG10 [53], were reported with oncogenic functions. Therefore, AC073103.1 might act as oncogene via competing for miR-150-5p.

Conclusions
Here we found the correlation of miR-150-5p and tumor in ltrating B cells and showed that miR-150-5p may improve the immune response in TME. This study provides insight into the cross-talk of HNSC cells and B cells in the TME through exosomal AC073130.1 and miR-150-5p.

Declarations
Ethics approval and consent to participate The study was approved by the institutional ethics committee of Sun Yat-sen University (Guangzhou, China).

Consent for publication
Not applicable.

Availability of data and materials
The information of this study here is obtained from GEO (https://www.ncbi.nlm.nih.gov/geo/) and the TCGA (https://portal.gdc.cancer.gov/).

Competing interests
The authors declare that they have no con icts of interest with the contents of this article.   Gene expression pro le associated with miR-150-5p in HNSC. A, miR-150-5p was down-regulated in HNSC. B-D, miR-150-5p expression was associated with overall survival and progress-free survival. E-G, Gene expression pro le associated with miR-150-5p. H-I, GO and KEGG analysis of different genes between miR-150-5p high expression group and low expression group in HNSC. J-N,GSEA analysis of different genes between miR-150-5p high expression group and low expression group in HNSC.