miR-873 and miR-105-2 May Affect the Tumour Microenvironment and are Potential Biomarkers for Lung Adenocarcinoma

Lung cancer, one of the most common malignant tumours, is the leading cause of cancer-related deaths worldwide and causes more than 1.5 million deaths every year.¹ Lung adenocarcinoma (LUAD), one of the predominant subtypes of non-small-cell lung cancer (NSCL), accounts for approximately 40% of all lung cancers.² It usually arises at more distal airways and is now viewed as the most heterogeneous and invasive subtype of lung cancer.³ Although the identification of oncogenic driver mutants and the clinical application of immunotherapy have improved LUAD treatment in recent years, the 5-year survival rate of LUAD is still less than 20%.⁴ Therefore, the identification of novel biomarkers is extremely important for the early diagnosis and treatment of LUAD.

MicroRNAs are endogenous small noncoding RNAs approximately 20nt long. They are able to regulate the expression of target genes by binding with their 3’ untranslated region (3ʹUTR).⁵ Although several microRNAs elevate the expression of target genes,⁶ the vast majority of microRNAs repress the expression of target genes by promoting mRNA degradation⁷ or inhibiting the translation process.⁸ By inhibiting the expression of tumour suppressor genes or oncogenes, microRNAs are able to promote or suppress the metastasis and prognosis of lung cancer.⁹ Moreover, microRNAs play an important role in the TME by regulating the growth and proliferation of immune cells and stromal cells.¹⁰

The ESTIMATE algorithm, developed by Yoshihara et al, estimates the relative ratios of stromal and immune components in the tumour microenvironment (TME) by calculating the StromalScore and ImmuneScore.¹¹ The CIBERSORT algorithm is able to calculate the abundance profiles of 22 kinds of immune cells in the TME.¹² Mining The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), Oncomine and other public databases by these two algorithms, BTK,¹³ TTC21A,¹⁴ JAK1,¹⁵ ADAM12 and ERG¹⁶ were found to be related to the TME and affected the prognosis of LUAD. These studies revealed that investigating the TME by ESTIMATE and the CIBERSORT algorithm was an effective way to screen out potential prognostic markers and therapeutic targets. However, most of these studies focused on protein-coding genes; there are few reports about the role of microRNAs in the TME of LUAD. Due to the important roles of microRNAs in the TME and LUAD, our study explored microRNAs related to the TME and the prognosis of LUAD by the ESTIMATE and CIBERSORT algorithms.

Materials and Methods

Data Source

Level 3 RNA-seq and microRNA-seq data of 570 LUAD samples, including 57 normal samples and 513 tumour samples, were downloaded from the TCGA database (https://portal.gdc.cancer.gov/, accessed June 30, 2021) along with their corresponding clinical data. GSE102286, containing expression information of microRNAs in 88 pairs of NSCLC and paracancerous tissues, was downloaded from GEO (https://www.ncbi.nlm.nih.gov/geo/). This study was approved by the ethics committee of Deyang People’s Hospital.

Calculation of the StromalScore and ImmuneScore

The StromalScore and ImmuneScore were positively correlated with the relative proportions of stromal and immune components in the TME. The ESTIMATEScore, the sum of the StromalScore and ImmuneScore, comprehensively reflects the relative proportion of stromal and immune components in the TME. The StromalScore, ImmuneScore and ESTIMATEScore were calculated by the ESTIMATE package of R language in our study.

Detection of Differentially Expressed microRNAs (DEMs) and Differentially Expressed Genes (DEGs)

LUAD samples were divided into a high-score group and a low-score group according to the median values of the StromalScore and ImmuneScore, respectively. Differential expression analysis between the high- and low-score groups was performed with the edgeR package.¹⁷ P < 0.05 and |log2FC| > 1 were used as the criteria to screen out DEMs and DEGs. Heatmaps were generated using the ggplot package of R language.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analyses

GO and KEGG enrichment analyses of DEGs were performed by the database for annotation, visualization, and integrated discovery (DAVID) (https://david.ncifcrf.gov/home.jsp).¹⁸ P < 0.05 was considered statistically significant.

Calculating the Abundance Profile of Tumour-Infiltrating Immune Cells (TICs) in the TME

The CIBERSORT algorithm, which analyses the relative expression levels of genes in each sample according to given gene expression profiles, was applied to estimate the TIC abundance profile in the TME.¹² Then, 417 samples with P < 0.05 were selected for subsequent analysis.

Prediction of microRNA Target Genes

Prediction of microRNA targets was performed by miRANDA (http://www.microrna.org/microrna/home.do),¹⁹ TargetScan (http://www.targetscan.org/vert_72/)²⁰ and Diana-microT-CDS (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=microT_CDS/index).²¹ Because miR-873 and miR-105-2 were downregulated in high StromalScore and ImmuneScore patients, upregulated DEGs, which were also predicted to be suppressed by miR-873 and miR-105-2, were screened out as potential targets of miR-873 and miR-105-2.

Statistical Analyses

R language and SPSS 19.0 (IBM, New York, USA) were used for statistical analysis. Differential analysis of transcriptomics data was performed using R language. Qualitative data were analysed by the chi-square test. Difference analyses of normally distributed and nonnormally distributed quantitative data were analysed by t tests and rank-sum tests, respectively. The Pearson and Spearman correlation tests were performed for correlation analyses of normally distributed and nonnormally distributed data, respectively. Kaplan–Meier survival analysis was performed to compare the survival times of different groups. P < 0.05 indicated that the difference was statistically significant.

Results

Relationship Between the StromalScore, ImmuneScore and Clinical Characteristics

A cohort of 513 LUAD patients from TCGA database were included in this study. The StromalScore and ImmuneScore of each patient were calculated by the ESTIMATE algorithm. The relationship between these scores and clinical characteristics was further analysed. Our results revealed that the StromalScore was related to the age (P = 0.001), sex (P = 0.009) and M classification (P = 0.007) of LUAD patients. Meanwhile, the ImmuneScore was associated with the age (P = 0.003), sex (P = 0.007), pathological stage (P = 0.039) and T classification (P = 0.003) of LUAD patients (Table 1).

Table 1 The Relationship Between StromalScore, ImmuneScore and Clinical Characteristics of LUAD Patients

Impact of the StromalScore, ImmuneScore and ESTIMATE Score on the Prognosis of LUAD

Kaplan–Meier survival analysis was performed to explore the impact of the StromalScore, ImmuneScore and ESTIMATEScore on the prognosis of LUAD. Then, a low StromalScore, ImmuneScore and ESTIMATEScore were found to decrease the overall survival time of LUAD patients (Figure 1A–C). These results revealed that the relative proportions of stromal and immune components affected the prognosis of LUAD. Therefore, screening out DEMs and DEGs between the high- and low-score groups was an effective way to identify biomarkers for LUAD.

Figure 1 The influence of StromalScore, ImmuneScore and ESTIMATEScore on the prognosis of LUAD. Kaplan–Meier survival curve of LUAD patients with high or low StromalScore (A) ImmuneScore (B) and ESTIMATEScore (C) Green and blue lines represented high- and low- score groups, respectively.

DEMs and DEGs Between the High- and Low-Score Groups

DEMs between the different StromalScore and ImmuneScore groups were screened out. There were 74 microRNAs differentially expressed between the low and high StromalScore groups. Among these microRNAs, 72 and 2 microRNAs were downregulated and upregulated in the high StromalScore group, respectively (Figure 2A and C). Meanwhile, 78 DEMs, including 67 downregulated and 11 upregulated microRNAs in high ImmuneScore patients, were screened out (Figure 2B and D). Taken together, 62 microRNAs were downregulated in both high StromalScore and high ImmuneScore patients (Figure 2E).

Figure 2 DEMs between high and low score groups. (A) Heatmap of DEMs between high and low StromalScore groups. Red and green represented upregulated and downregulated microRNAs, respectively. (B) Heatmap of DEMs between high and low ImmuneScore groups. Red and green represented upregulated and downregulated microRNAs, respectively. (C) Volcano plot of DEMs between high and low StromalScore groups. Red and green dots represented up-regulated and down-regulated DEMs, respectively. Black dots represented microRNAs without significantly changed. (D) Volcano plot of DEMs between high and low ImmuneScore groups. Red and green dots represented up-regulated and down-regulated DEMs, respectively. Black dots represented microRNAs without significantly changed. (E) Venn plot of down-regulated DEMs shared by StromalScore and ImmuneScore.

DEGs were also screened out between different StromalScore and ImmuneScore groups. The expression levels of 1947 and 1436 genes were increased and decreased in high StromalScore patients, respectively (Figure 3A and C). In high ImmuneScore patients, the expression of 1509 and 1684 genes was increased and decreased, respectively (Figure 3B and D). Taken together, 669 and 946 DEGs were upregulated and downregulated in both high StromalScore and high ImmuneScore patients, respectively (Figure 3E and F). GO enrichment analysis revealed that these DEGs were mainly involved in biological processes including immune response, regulation of immune response and adaptive immune response (Figure 3G), while KEGG enrichment analysis revealed that these DEGs were mainly involved in signal pathways such as systemic lupus erythematosus, neuroactive ligand–receptor interaction, cytokine–cytokine receptor interaction and alcoholism (Figure 3H).

Figure 3 DEGs between high and low score groups. (A) Heatmap of DEGs between high and low StromalScore groups. Red and green represented genes with high and low expression levels, respectively. (B) Heatmap of DEGs between high and low ImmuneScore groups. Red and green represented genes with high and low expression levels, respectively. (C) Volcano plot of DEGs between high and low StromalScore groups. Red and green dots represented up-regulated and down-regulated DEGs, respectively. Black dots represented genes without significantly changed. (D) Volcano plot of DEGs between high and low ImmuneScore groups. Red and green dots represented up-regulated and down-regulated DEGs, respectively. Black dots represented genes without significantly changed. (E) Venn plot of up-regulated DEGs shared by StromalScore and ImmuneScore. (F) Venn plot of down-regulated DEGs shared by StromalScore and ImmuneScore. (G) GO enrichment result of DEGs shared by StromalScore and ImmuneScore. The size and color of each bubble represented gene number and p value of each term, respectively. (H) KEGG enrichment result of DEGs shared by StromalScore and ImmuneScore. The size and color of each bubble represented gene number and p value of each term, respectively. Abbreviations: BP, CC and MF, biological process, cellular component and molecular function, respectively.

miR-873 and miR-105-2 Were Associated with the Prognosis and Occurrence of LUAD

Survival analysis was performed to explore the effect of DEMs on the prognosis of LUAD. Increased expression levels of miR-873, miR-105-2 and miR-516a-2 significantly shortened the survival time of LUAD patients (Figure 4A–C). A paired rank-sum test was performed to investigate the expression levels of these microRNAs in paracancerous and LUAD tissues. Among these 3 DEMs, only the expression levels of miR-873 and miR-105-2 increased significantly in tumour tissues (Figure 4D and E). Therefore, we verified the expression patterns of miR-873 and miR-105-2 using GSE102286. Consistent with the aforementioned results, the expression levels of miR-873 and miR-105 increased significantly in NSCL tissues (Figure 4F and G). Taken together, these results revealed that miR-873 and miR-105-2 were closely related to the prognosis and occurrence of LUAD.

Figure 4 Survival analysis, difference analysis and functional influence of miR-873 and miR-105-2 on LUAD. Kaplan–Meier survival curve of LUAD patients with high or low expression level of miR-873 (A), miR-105-2 (B) and miR-516a-2 (C). Green and blue lines represent patients with high- and low- expression levels, respectively. Paired rank-sum analysis results of the expression levels of miR-873 (D) and miR-105-2 (E) in LUAD and paracancerous tissues. Red and blue rectangles represented cancerous and paired paracancerous tissues, respectively. Paired rank-sum analysis results of the expression levels of miR-873 (F) and miR-105 (G) in NSCLC and paracancerous tissues. Red and blue rectangles represented cancerous and paired paracancerous tissues, respectively. (H) Predicted targets of miR-873 and miR-105-2. Red rectangles and green triangles represented DEGs upregulated in high score groups and DEMs downregulated in high score groups, respectively. (I) KEGG enrichment result of potential targets of miR-873 and miR-105-2. The size and color of each bubble represented gene number and p value of each term, respectively. (J) GO enrichment result of potential targets of miR-873 and miR-105-2. The size and color of each bubble represented gene number and p value of each term, respectively. Abbreviations: BP, CC and MF, biological process, cellular component and molecular function, respectively.

The potential targets of miR-873 and miR-105-2 were predicted to explore their functional effect. Then, 109 and 126 DEGs were predicted to be regulated by miR-873 and miR-105-2, respectively (Figure 4H). KEGG analysis revealed that these target genes were involved in cytokine–cytokine receptor interactions, haematopoietic cell lineages, intestinal immune networks for IgA production, and systemic lupus erythematosus (Figure 4I). Meanwhile, GO enrichment analysis results revealed that these target genes were mainly enriched in immune response, adaptive immune response and inflammatory response (Figure 4J). Interestingly, 17 out of 28 KEGG pathways modulated by miR-873 and miR-105-2, including cytokine–cytokine receptor interaction, haematopoietic cell lineage, intestinal immune network for IgA production and systemic lupus erythematosus, were related to DEGs between different StromalScore and ImmuneScore groups. These results suggested that miR-873 and miR-105 may affect the TME via these pathways.

miR-873 and miR-105-2 May Alter the Relative Proportions of TICs

The CIBERSORT algorithm and correlation analysis were used to analyse the relative proportions and the relationships between 22 kinds of TICs in LUAD patients (Figure 5A and B). Rank-sum tests were performed to compared the relative proportions of 22 kinds of TICs in cancerous and paracancerous tissues. Our results revealed that the relative proportions of most TICs in cancerous and paracancerous tissues were significantly different, implying the important roles of TICs in LUAD (Figure 5C). Then, correlation and difference analyses were performed to explore the effects of miR-873 and miR-105 on TICs. Follicular helper T cells, resting NK cells, neutrophils and CD4 naive T cells, monocytes, resting dendritic cells, and resting mast cells were positively and negatively correlated with the expression of miR-873, respectively (Figure 6A). Meanwhile, high expression of miR-873 elevated and decreased the relative proportions of follicular helper T cells, neutrophils and memory resting CD4 T cells, monocytes, and resting dendritic cells, respectively (Figure 6B). Subsequent analysis revealed that miR-105-2 was positively and negatively correlated with the abundance of resting NK cells, M0 macrophages, M1 macrophages and memory resting CD4 T cells, monocytes, M2 macrophages, resting dendritic cells, activated dendritic cells, and resting mast cells, respectively (Figure 6C). Additionally, a high expression level of miR-105-2 decreased and elevated the relative proportions of resting memory CD4 T cells, monocytes, resting dendritic cells and activated memory CD4 T cells, M0 macrophages, and M1 macrophages, respectively (Figure 6D). Intersection analysis results revealed that elevated expression levels of miR-873 and miR-105-2 decreased the relative proportions of monocytes and resting dendritic cells in the TME (Figure 6E). Taken together, our results implied that miR-873 and miR-105-2 may affect the immune activity of the TME by altering the relative proportions of TICs.

Figure 5 The abundance profile of TICs in LUAD samples and correlation analysis. (A) Bar plot showing the proportion of 22 kinds of TICs in LUAD samples. Each column represented a LUAD sample. Different colors represented different kinds of TICs. (B) Heatmap showing the correlation between 22 kinds of TICs. The color of each box represented the correlation index between two kinds of TICs. (C) Bar plot showing the proportion of 22 kinds of TICs in cancerous and paracancerous tissues. NS, *, ** and *** represented P>0.05, P<0.05, P<0.01 and P<0.001, respectively. Red and blue rectangles represented cancerous and normal tissues, respectively.

Figure 6 Correlation of proportions of TICs and the expression levels of miR-873 and miR-105-2. (A) Scatter plots indicated the correlation between the expression level of miR-873 and 7 kinds of TICs. The blue line in each plot was fitted linear model indicating the correlation between indicated TIC and the expression level of miR-873. (B) Box plot indicated the influence of high and low expression level of miR-873 on proportions of 22 kinds of TICs. NS and * represented P>0.05 and P<0.05, respectively. (C) Scatter plots indicated the correlation between the expression level of miR-105-2 and 9 kinds of TICs. The blue line in each plot was fitted linear model indicating the correlation between indicated TIC and the expression level of miR-105-2. (D) Box plot indicated the influence of high and low expression level of miR-105-2 on proportions of 22 kinds of TICs. NS, * and ** represented P>0.05, P<0.05 and P<0.01, respectively. (E) Venn plot of TICs correlated with the expression levels of miR-873 and miR-105-2 codetermined by difference and correlation tests.

Discussion

The TME is composed of immune cells, stromal cells, endothelial cells, extracellular matrix, blood vessels and other kinds of cells. Although immunotherapy using immune checkpoint inhibitors such as antibodies against PD-1 and CTLA-4 has improved the prognosis of NSCL in recent years, the efficacy of immunotherapy is affected by the TME.²² Moreover, accumulating evidence has revealed that the TME plays an important role in the occurrence, metastasis, recurrence and treatment of tumours.²³ Therefore, investigating TME-related genes will help to further reveal the molecular mechanism of LUAD and identify potential biomarkers.

In our study, the effect of the TME on the prognosis of LUAD was investigated by mining the TCGA database. Consistent with previous reports,^13,24,25 our results showed that high values for the StromalScore, ImmuneScore and ESTIMATEScore improved the prognosis of LUAD. However, most of these researches focused on protein-coding genes. Although microRNAs played an important role in the TME, few research explored their function in the TME. Therefore, we explored the effect of microRNAs in the TME of LUAD via bioinformatics analysis. In both the high StromalScore and high ImmuneScore groups, 62 DEMs were screened out to be significantly downregulated. Among these microRNAs, miR-873, miR-105-2 and miR-516a-2 were associated with poor prognosis of LUAD. Moreover, the expression levels of these microRNAs in cancerous and paired paracancerous tissues were verified using a TCGA LUAD dataset and an independent GEO NSCL dataset. Compared with paracancerous tissues, the expression levels of miR-873 and miR-105-2 increased significantly in cancerous tissues. Finally, our bioinformatics analysis results revealed that miR-873 and miR-105-2 may alter the abundance profiles of TICs in the TME and then affect the occurrence and prognosis of LUAD.

MiR-873, located on chromosome 9p21.1, was able to affect the proliferation, invasion and metastasis of many cancers by inhibiting the expression of downstream targets.²⁶ Recent studies revealed that miR-873 was upregulated in LUAD and promoted the proliferation, migration and progression of tumour cells.^27,28 Additionally, elevated expression of miR-873 significantly shortened the survival time of LUAD patients.²⁹ MiR-105-2, residing on chromosome Xq28, belongs to the miR-767-105 cluster and leads to poor prognosis of hepatocellular carcinoma.^30,31 A high plasma miR-105 level resulted in poor prognosis of NSCL and was a potential biomarker for early diagnosis of NSCL.³² In addition, miR-105 was inhibited by LINC00261 and affected the metastasis and proliferation of NSCLC.³³ Consistent with previous reports, our results revealed that miR-873 and miR-105-2 shorten the survival time of LUAD patients by altering the relative proportions of TICs in the TME.

Growing evidence suggests that TICs affect tumour progression and the immunotherapeutic response.³⁴ Therefore, we explored the influence of miR-873 and miR-105-2 on TICs in our study. Our results revealed that miR-873 and miR-105-2 were correlated with the relative proportions of several kinds of TICs. Among these kinds of TICs, we focused on monocytes and resting dendritic cells, which were more closely associated with the expression levels of miR-873 and miR-105-2. Monocytes are a heterogenic population of mononuclear phagocytes that comprise a major population of innate immune cells. By inducing tumour cytotoxicity, suppressing metastasis, engulfing tumour materials and negatively regulating Treg cells, monocytes were able to suppress tumour progression.³⁵ Although the proportions of CD14+ monocytes were similar in tumour and normal tissues, the proportion of CD16+ monocytes decreased significantly in NSCL tumour tissues.³⁶ Dendritic cells are a group of antigen-presenting cells involved in adaptive immune responses. Resting dendritic cells, also known as immature dendritic cells, are primarily localized in peripheral tissues and can capture antigens efficiently. By regulating T-cell activation, resting dendritic cells are mainly involved in maintaining self-tolerance.³⁷ Previous research revealed that immature dendritic cells led to a better prognosis of lung cancer.³⁸ Later studies unveiled that dendritic cells were more likely to remain immature in the TME in NSCLC.^39,40 These results implied that monocytes and resting dendritic cells may suppress the progression of lung cancer. Our results revealed that miR-873 and miR-105-2 were negatively correlated with the abundance of monocytes and resting dendritic cells in the TME of LUAD. Taken together, our results implied that miR-873 and miR-105-2 may suppress the progression of NSCLC by decreasing the relative proportions of monocytes and resting dendritic cells.

Conclusion

In conclusion, our results revealed that miR-873 and miR-105-2 were closely related to the TME and may affect the prognosis of LUAD by altering the relative proportions of TICs. Although these findings still need to be further investigated in vitro and in vivo, our results identified 2 potential biomarkers for LUAD.