Roles of the m6A Modification of RNA in the Glioblastoma Microenvironment as Revealed by Single-Cell Analyses.

Purpose: Glioblastoma multiforme (GBM) is a common and aggressive form of brain tumor. The N6-methyladenosine (m6A) mRNA modification plays multiple roles in many biological processes and disease states. However, the relationship between m6A modifications and the tumor microenvironment in GBM remains unclear, especially at the single-cell level. Experimental Design: Single-cell and bulk RNA-sequencing data were acquired from the GEO and TCGA databases, respectively. We used bioinformatics and statistical tools to analyze associations between m6A regulators and multiple factors.
Results: HNRNPA2B1 and HNRNPC were extensively expressed in the GBM microenvironment. m6A regulators promoted the stemness state in GBM cancer cells. Immune-related BP terms were enriched in modules of m6A-related genes. Cell communication analysis identified genes in the GALECTIN signaling network in GBM samples, and expression of these genes (LGALS9, CD44, CD45, and HAVCR2) correlated with that of m6A regulators. Validation experiments revealed that MDK in MK signaling network promoted migration and immunosuppressive polarization of macrophage. Expression of m6A regulators correlated with ICPs in GBM cancer cells, M2 macrophages and T/NK cells. Bulk RNA-seq analysis identified two expression patterns (low m6A/high ICP and high m6A/low ICP) with different predicted immune infiltration and responses to ICP inhibitors. A predictive nomogram model to distinguish these 2 clusters was constructed and validated with excellent performance.
Conclusion: At the single-cell level, m6A modification facilitates the stemness state in GBM cancer cells and promotes an immunosuppressive microenvironment through ICPs and the GALECTIN signaling pathway network. And we also identified two m6A-ICP expression patterns. These findings could lead to novel treatment strategies for GBM patients.

Introduction

Glioblastoma multiforme (GBM) is the most common primary tumor of the central nervous system, and it is extremely aggressive (1). Standard treatment for GBM is surgical resection followed by chemoradiotherapy, with a 5-year survival of 7.2% (1, 2). Immunotherapy has achieved growing success across systemic cancers, and has become a prominent player in the treatment of GBM (3). Nonetheless, the high degree of genetic heterogeneity and immunosuppressive microenvironment that characterize GBM represent important challenges to the application of immunotherapy to this disease (1). Advancements in knowledge of the immune cells in the GBM microenvironment, particularly glioma-associated microglia, macrophages and T cells, might lead to novel strategies to strengthen anti-tumor immunity (3).

An important aspect of the behavior of these key immune cells involves the N6-methyladenosine (m6A) modification, which has emerged as the most abundant chemical modification of protein-coding and noncoding RNAs (4, 5). To date, regulators of the m6A modification have been reported to be involved in cancer biology, including cancer progression and other processes (4). Notably, the roles of m6A regulators in mediating immunotherapy resistance have been highlighted in recent studies. For instance, expression of METTL3 in macrophages was suggested to synergize with PD-1-based therapy in B16 melanoma (6), and METTL3 and METTL14 have been shown to regulate immune response to anti-PD-1 treatment in melanoma and colorectal carcinoma (7). Together, these findings indicate a role for m6A regulators as potential therapeutic targets in anticancer immunotherapy.

Several m6A regulators have been reported to be upregulated and to play vital roles in GBM, including enhancing cell self-renewal and proliferation in glioblastoma stem cells (GSCs) (8). However, some opposing findings have shown lower expression levels of some m6A writers in GBM and have detected anticancer properties of m6A regulators (8, 9). This conflict suggests a complex role of m6A-related methylation in the occurrence and development of GBM. Notably, though, FTO, YTHDF2, and RBM15 were found to have prognosis predictive value in GBM (8, 10). In addition, patients with higher m6A related risk scores were more sensitive to temozolomide treatment and showed lower drug resistance overall (10).

These findings suggested critical and complex roles of m6A modification in GBM, and more research opportunities and challenges have emerged. Recent development of single-cell analysis methods provides more comprehensive approaches to explore the potential mechanism of action of m6A modifications in GBM at the cellular level. Herein, for the first time, we used scRNA-seq data to thoroughly analyze () the roles of m6A modifications in the GBM microenvironment, especially its relationship with functional states, cell communication and immune checkpoints (ICPs).

Figure 1

Schematic diagram of the study design.

Materials and Methods

Datasets

Five publicly available human GBM scRNA-seq datasets (11–15) were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Bulk RNA-seq datasets and DNA methylation profiles were obtained from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and the Repository of Molecular Brain Neoplasia Data (REMBRANDT) (16) datasets. These bulk datasets were downloaded via the University of California, Santa Cruz (UCSC) Xena browser (https://xenabrowser.net/) and Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/) websites. A summary description of included scRNA-seq datasets is provided in .

Single-Cell RNA-seq Analysis

The “Seurat” R package (version 4.0.2) (17) was used to perform single-cell RNA-seq analyses. We only included in this study samples with at least 10 000 detected genes. Quality control (QC) was based on following standards: 1) genes detected in fewer than 3 cells were excluded; 2) cells with fewer than 200 total detected genes were excluded; 3) cells with at least 5% (for GSE141383, GSE84465, GSE103224, GSE89567) or 10% (for GSE138794) of mitochondrial genes were excluded; 4) cells with at most 3% (for GSE141383, GSE138794, GSE103224, GSE89567) or 1% (for GSE84465) of ribosomal genes were excluded; 5) cells with at least 0.1% of hemoglobin genes were excluded and 6) cells with total gene expression between first quartile − 1.5 × interquartile range (IQR) and third quartile + 1.5 × IQR were retained to exclude cellular doublets.

Unnormalized datasets were normalized with the “NormalizeData” function in R. Then, these five datasets were merged together, and the zero imputation method was applied to investigate missing data relating to significant markers (). The batch effect was removed via the “fastMNN” algorithm from the “SeuratWrappers” R package (version 0.3.0).

Principal component analysis (PCA) was performed, and the top 20 PCs were used for cluster classification. The initial 20 PCs were also utilized in the further visualization process through the t-distributed stochastic neighbor embedding (t-SNE) algorithm. After cluster classification, different cell clusters were identified and annotated manually within the CellMarker database (http://biocc.hrbmu.edu.cn/CellMarker/) (18). The cell cycle score was calculated with the “CellCycleScoring” function. Cell communication network analysis was conducted with the “CellChat” R package (version 1.1.1). Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) was used to compare the expression levels of m6A regulators between GBM and normal samples with bulk RNA-seq datasets.

Single-Sample Gene Set Enrichment Analysis

The ssGSEA algorithm was conducted with the “GSVA” R package (version 1.34.0) to calculate enrichment scores for several gene signatures. The gene set for m6A included 13 m6A readers (YTHDC1, YTHDC2, IGF2BP1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF2, YTHDF3, HNRNPA2B1, HNRNPC, ELAVL1, LRPPRC and FMR1), 8 m6A writers (METTL3, METTL14, WTAP, KIAA1429, RBM15, RBM15B, ZC3H13 and CBLL1) and 2 m6A erasers (FTO and ALKBH5). The gene lists () for other biological states were from CancerSEA (http://biocc.hrbmu.edu.cn/CancerSEA).

m6A Related Genes

Analyses of Pearson correlations between m6A scores and gene expression were conducted to filter m6A related genes under the following criteria: (1) R value greater than 0.3 or less than -0.3 and (2) P < 0.05. In M2 macrophages, we strengthened the criteria for m6A related genes to R value greater than 0.5 or less than -0.5. The Molecular Complex Detection (MCODE) plugin from Cytoscape software (version 3.8.1) was applied to identify core gene modules with scores greater than 10.

Functional Annotation

Enrichment analyses of biological processes (BP) Gene Ontology (GO) annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed based on the “clusterProfiler” R package (version 3.14.3) with adjusted P < 0.05 as the cutoff criterion.

Protein-Protein Interaction Network Analysis

The PPI network was constructed with Search Tool for the Retrieval of Interacting Genes (STRING) (https://string-db.org/) and visualized with Cytoscape.

Bulk RNA-seq Analysis

To better understand the interaction between m6A regulators and immune checkpoints (ICPs), the “ConsensusClusterPlus” R package (version 1.50.0) was applied to identify distinct clusters of bulk RNA-seq samples with correlated expression levels of m6A regulators and ICPs. The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was utilized to predict potential therapeutic responses of ICP inhibitors. Morpheus (https://software.broadinstitute.org/morpheus) was used to identify differentially expressed genes (DEGs) between clusters. The cutoff criteria for DEGs were: P < 0.05 and |log2 fold-change| ≥ 1. An enrichment analysis of DEGs was conducted and hub genes in DEGs were visualized in a PPI network. The ssGSEA for various gene signatures was performed for each sample. We also conducted the Cell Type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) analysis of the abundance scores of immune cells in GBM samples. A Pearson correlation analysis among the methylation levels and expression levels of m6A regulators and ICPs was performed and visualized with Cytoscape software.

Construction of Predictive Model

To sort GBM patients into 2 clusters, the TCGA and REMBRANDT datasets were used as a training dataset and validation dataset separately. A least absolute shrinkage and selection operator (LASSO) regression analysis was used to identify all genes with P < 0.05. Then, best subset selection was performed to determine the final multivariate logistic regression model. A nomogram was constructed to show the predictive model. The evaluation of the nomogram model was performed with receiver operating characteristic (ROC) curve analysis and calibration curve analysis. “glmnet” (version 4.1), “rms” (version 6.1.0), and “timeROC” (version 0.4) R packages were used for the construction of predictive model.

Sample Source

Tumor tissues were surgically harvested at Jinling Hospital in accordance with institution‐approved protocols. The tissues were collected from 2 glioma patients (without preoperative radiotherapy or chemotherapy), confirmed pathologically as specimens of grade IV gliomas, and fixed in formalin before paraffin-embedded. Tissues from two nontumoral cases were also collected as control. Written informed consent had been obtained from participants.

Cell Lines and Culture

The human GBM cell line U87MG were purchased from cell bank of Shanghai Institute of Biochemistry and Cell Biology and grown in Dulbecco’s Modified Eagle’s medium (DMEM) containing 10% fetal bovine serum (FBS, Gibco, USA). THP-1 was kindly provided by Wentao Liu laboratory, Nanjing Medical University and cultured in 1640 medium with 10% FBS. Cells were cultured at 37°C in a humidified atmosphere with 5% CO2.

Proteins and Reagents

Antibodies against MDK (11009-1-AP) and CD206 (60143-1-lg) was purchased from proteintech, and that against LRP1 (ab92544) was from Abcam. The β-Actin polyclonal antibody (AP0060) and goat anti-rabbit IgG-HRP (BS13278) were purchased from Bioworld. Midkine protein (abs00930) was from Absin Bioscience Inc. The jetPRIME transfection reagent (101000046) was from Polyplus. Anti-CD206 APC (17–2069–41) and its isotype control (17–4714–81) were purchased from eBioscience. And anti-CD11b FITC (101205), its isotype control (400633), and cell staining buffer (420201) was purchased from BioLegend.

Transfections With siRNA

The MDK small interfering RNA (siMDK; forward, 5’-GACCA AAGCA AAGGC CAAATT-3’; reverse, 5’-UUUGG CCUUU GCUUU GGUCTT-3’) and negative control siRNA (siCon; forward, 5’-UUCUC CGAAC GUGUC ACGUTT-3’; reverse, 5’-ACGUG ACACG UUCGG AGAATT-3’) were synthesized by GenePharma (Shanghai, China). U87MG was transfected with 100ul/L MDK siRNA and negative control siRNA respectively, using jetPRIME transfection reagent. Cells were transfected with siRNAs for 48 h prior to the following experiments.

Western Blot Assay

Total proteins of cells were extracted, quantified using BCA kit (23235, Thermo), separated by 10% SDS-PAGE, and transferred onto polyvinylidene difluoride membrane. The membrane was blocked with 5% defatted milk powder for 1.5 h, and was incubated overnight in primary antibody (anti-β-Actin at 1:10000; anti-MDK at 1:1000) at 4°C. On the following day, the membrane was incubated in secondary antibody (goat anti-rabbit IgG-HRP 1:20000) for 1 h at room temperature. After adding ECL chromogenic substrate (Millipore, US), the membrane was imaged using a gel imaging system.

qPCR

RNA was extracted using Ultrapure RNA Kit (CW0597, CWBIO). cDNA was synthesized using a reverse transcription kit (AE311, TransGern) following the instructions in the manual. The endogenous levels of MDK mRNA were determined using the SYBRGreen PCR Kit (AQ131, TransGern). PCR conditions were: denaturation at 94°C for 30 s, followed by 40 cycles of denaturation at 94°C for 5 s, annealing at 61°C for 35 s, followed by elongation at 95°C for 10 s. GAPDH was used as reference gene for MDK. The following primers were used: MDK-F (5’-AAGGATTGCGGCGTGGGTTTC-3’), MDK-R (5’-TGGCGGACTTTGGTGCCTGTG-3’), GAPDH-F (5’-AATCCCATCACCATCTTCCA-3’), and GAPDH-R (5’-AAATGAGCCCCAGCCTTCT-3’). Relative expression level of MDK was calculated using 2-▵▵Ct approach.

Transwell Assay for Cell Migration

THP-1 was primed with 185 ng/ml PMA for 24 hours, and then added to the top chamber of Transwell unit at a density of 3 × 105 cells/unit. Next, U87MG cell suspension was added into the lower chamber. The transwell chamber was cultured at 37°C for 24 h in a cell incubator. The cells were fixed with methanol for 15 min, and stained with 0.1% crystal violet solution for 5 min. Then, cells on the upper surface of the chamber bottom were removed gently with a cotton swab. Finally, observe the invaded cells on the lower surface under an inverted microscope. Three separate membranes were analyzed for each condition.

Flow Cytometry

Staining for cell surface markers was done by re-suspending each sample in 100 µl cell staining buffer containing the antibody cocktail. Three groups were set here. Isotype controls were added into group 1. Anti-CD206 APC and anti-CD11b FITC were added into group 2 (siCon) and 3 (siMDK). Cells were incubated at 4°C for 30 minutes and then washed with cell staining buffer. Cells with CD11b+/CD206+ phenotype were identified as M2 macrophages. Data were immediately acquired using CytoFLEX (Beckman Coulter) and analyzed with FlowJo software (version 10.5.3).

Immunohistochemistry

Paraffin‐embedded clinical tissue specimens were sectioned, dewaxed, and dehydrated. Antigen retrieval was conducted using a pressure cooker for 3 min. Then, sections were washed with 3% methanol H2O2. Subsequently, the sections were incubated overnight at 4°C using primary antibodies against MDK (11009-1-AP, 1:500) and LRP1 (ab92544, 1:300), and treated with biotin-free EnVision detection kit (Dako) for secondary antibody. The labeled antigens were visualized by 3,3’-diaminobenzidine tetrahydrochloride as a chromogen. Finally, the sections were counterstained with hematoxylin. ImageJ software and IHC Profiler plugin were applied for qualitative assessment of IHC slides. Five random fields were observed under microscope.

Sections were scored using Histochemistry score (H-SCORE) method. Staining proportion was scored from 0 to 3, with 0 if negative, 1 if < 25% of cells stained positively, 2 if 26–50% of cells stained positively, and 3 if > 50% of cells stained positively. In addition, the staining intensity was scored as 0, no staining; 1, weak; 2, moderate; and 3, strong. These two scores were multiplied by each other to calculate the expression score of MDK and LRP1, with 0-3 represents low expression, and 4–9 represents high expression.

Immunofluorescence Staining

Segments of samples were dehydrated with 15% and 30% sucrose solution sequentially, embedded in OCT compound, and cryosectioned into 8 μm sections. Fixation was conducted with 4% paraformaldehyde for 15 min. After washing with PBS at room temperature, samples were incubated with 0.1% Triton X-100 for 5-15 min, were followed by washing with PBS. Sections were blocked with 10% goat serum for 1 h at room temperature. Primary antibody was added for overnight incubation at 4°C.

On the following day, primary antibody was removed by rinsing with PBS, and TRITC/FITC -labelled secondary antibody (1:100) was added and incubated at room temperature for 30 min and subsequently stained with DAPI (1:400) for 1-2 min. Coverslips were mounted. Images were obtained under a laser scanning confocal microscope.

Statistical Analysis

All statistical analyses were performed in R software (version 3.6.0 and 4.0.2). Student’s t-test and chi-square test or Fisher’s exact test were used for comparisons of continuous variables and categorical variables. Bar charts were drawn to display the comparisons using the “ggplot” R package (version 3.3.3). The Kaplan–Meier approach was performed to estimate survival, and the overall survival (OS) was compared with log-rank tests. We also applied Pearson correlation analyses to evaluate two continuous variables, after calculating the mean value of these variables in each sample. To avoid the influence of missing data, we excluded samples with missing data in correlation analyses, and only conducted the correlation analysis when there were more than 3 samples. The correlation analysis results were visualized with “corrgram” (version 1.13) and “ggplot” R packages. P < 0.05 values were set as indicating statistical significance.

Results

Analysis of scRNA-seq Data Identifies 7 Types of Cells in GBM

We obtained scRNA-seq profiles of 52 851 single cells from 22 GBM samples (). Following the application of a QC standard, 47 978 single cells were retained (; ). To visualize the distribution of the scRNA-seq profile, we employed a t-distributed stochastic neighbor embedding (t-SNE) algorithm to reduce the dimensionality of these datasets (). After appropriate elimination of the batch effect within these datasets with the “fastMNN” algorithm, the data was well integrated (). With unsupervised classification, cells were successfully classified into 17 clusters (). Based on the expression patterns of markers from CellMarker, we manually annotated these clusters as the following 7 cell types: 1) GBM cancer cells (expressing SOX2, PARP1, and CCND2); 2) M1 macrophages (CD68, CD74, TSPO, and CD86); 3) M2 macrophages (CD68, CD74, and CD163); 4) T cells or natural killer (NK) cells (CXCR4 and S100A4); 5) endothelial cells (A2M and APOLD1); 6) astrocytes (GFAP and SOX9); 7) oligodendrocytes (CNP, MBP, and PLP1) (). We also detected heterogeneous cell compositions among these included samples (). DEGs for cell types were identified (; ), and enrichment analyses for these DEGs were conducted to show related BP and pathways of each cell type ().

Figure 2

Identification of 17 cell clusters and 7 types of cells in GBM tumors. (A, B) tSNE plot of GBM cells before (A) and after (B) batch effect elimination. (C) Unsupervised classification successfully identified 17 cell clusters. (D) All 17 clusters were annotated by CellMarker according to the composition of the marker genes. (E) tSNE plot of 7 cell types. (F) Expression levels of marker genes in the 7 cell types. (G) Heatmap of differentially expressed features in each cell type. (H) Distribution of 7 cell types for all included cells. (I) Distribution of 7 cell types in each included sample.

HNRNPA2B1 and HNRNPC Show Extensive Expression in the GBM Microenvironment

First, we obtained the expression levels from the bulk transcriptome of 23 m6A regulators in GBM and normal samples from GEPIA, and found that most m6A regulators have relatively higher expression levels in GBM compared with that in normal samples. Especially notable were HNRNPA2B1, HNRNPC, and WTAP (). Then, to explore the expression patterns of m6A regulators in GBM samples at the single-cell level, we drew t-SNE plots of these m6A regulators. In general, all 23 m6A regulators were expressed in all 7 cell types, but not in all cells (; ).

Figure 3

Distribution of m6A regulators in GBM microenvironment. (A) Expression levels of m6A regulators in the 7 cell types. (B–O) tSNE plots of 14 m6A regulators. (P) tSNE plot of cells in 3 cell cycles. (Q) Distribution of cells in 3 stages of the cell cycle in 7 cell types. (R) Expression levels of m6A regulators in 3 stages of the cell cycle.

Notably, HNRNPA2B1 and HNRNPC showed extensive expression pattern in all cell types, but the expression of these m6A regulators was higher in M1 macrophages, GBM cancer cells, and T/NK cells. We further explored the expression levels of these m6A regulators in different cell cycle phases (). Heterogeneity of cell cycle distribution was observed in these included samples (). The GBM cancer cells exhibited the most and largest proportion of cells at S and G2M phases. Interestingly, other types of cells, except M1 macrophages, had about 30% to 50% cells in the S and G2M phases. In general, cells at G1 phase showed higher expression of m6A regulators than did cells at S and G2M phases (). For cells at S phase, HNRNPC had highest expression level relative to other m6A regulators. A higher expression level of HNRNPA2B1 was observed in cells at G2M phase. For cells at G1 phase, m6A regulators with higher expression levels were HNRNPA2B1, HNRNPC, WTAP, and YTHDC1.

Similar results were observed in all 7 types of cells, except for M2 macrophages (). For M2 macrophages at S phase, HNRNPA2B1 was more highly expressed than was HNRNPC. These results indicated intratumor heterogeneity of m6A regulators and intertumor heterogeneity of cell types and cell cycle phases in different samples. The following analyses of expression of m6A modification factors in GBM were based on subsets of the noted cell types.

m6A Regulators Associate With Functional States

To assess the potential relationship between m6A score and scores representing 16 functional states, we conducted a series of Pearson correlation analyses in all 7 cell types (). Specifically, we found that in GBM cancer cells, the m6A score was significantly correlated with the stemness score (R = 0.49, P < 0.05). To further explore the relationship between m6A regulators and stemness-related genes, we calculated the mean expression levels of genes in each sample, conducted Pearson correlation and PPI analyses between these two sets of genes (), and visualized these results into an m6A regulation network (). These results displayed extensive positive correlations between m6A regulators and stemness-related genes. Notably, SOX2, which can reprogram differentiated glioma cells to glioma stem-like cells (GSCs), has been identified as a bona fide m6A target of METTL3 ( 19), but the downstream m6A reader of SOX2 in GBM remains unknown. In colorectal carcinoma, IGF2BP2 has been proven to be a downstream m6A reader of SOX2 ( 20). Our regulation network suggested that both IGF2BP2 and ELAVL1 were potential downstream m6A readers of SOX2 in GBM ().

Figure 4

Relationship between m6A modification and 16 functional states. (A, B) Correlation heatmap (A) and correlation analysis (B) of m6A score and 16 functional state scores in 7 cell types. (C, D) Correlation heatmap (C) and correlation analysis (D) of m6A regulators and partial stemness related genes in GBM cancer cells. (E) Correlation network of m6A regulators and stemness related genes in GBM cancer cells. NS: P > 0.05, *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001.

Beside the stemness score in GBM cancer cells, we also found that differentiation score in GBM cancer cells; apoptosis, ferroptosis, and pyroptosis scores in M1 macrophages; epithelial-mesenchymal transition (EMT), invasion, and metastasis scores in M2 macrophages; apoptosis, DNA damage, and pyroptosis scores in T/NK cells; and angiogenesis, invasion, and pyroptosis scores in endothelial cells were all significantly correlated with m6A score. We further conducted Pearson correlation analyses between these gene sets and m6A regulators. By combining of these correlation analyses with PPI analyses, we constructed 13 m6A regulation networks (). These networks provided valuable information about potential regulatory mechanisms of m6A regulators in the GBM microenvironment, which will be explored and verified in further studies.

Modules of m6A Related Genes Associate With Immune-Related Terms

By conducting the Pearson correlation analysis between m6A score and all genes in all 7 cell types, we filtered m6A related genes for each cell type under the P < 0.05 criterion. With PPI analysis and the MCODE plugin from Cytoscape software (Degree cutoff = 2; Node score cutoff = 0.2; K-core = 2 and Max depth = 100), we identified functional modules of m6A related genes in a PPI network (; ). Enrichment analysis was applied for each cluster to explore potential related BP and pathways (). As expected, all cell types had clusters related to mRNA modification. Also, we found energy conversion and cell cycle related terms in clusters of most cell types. Specifically, terms related to antigen processing and presentation process were enriched in cluster 2 of GBM cancer cells and cluster 1 of endothelial cells (). In KEGG analyses, terms relating to “PD-L1 expression and PD-1 checkpoint pathway in cancer” were enriched in cluster 2 of M2 macrophages and cluster 1 and 2 of endothelial cells ().

m6A Regulators Promote Immunosuppressive Activities Through the GALECTIN Signaling Network

We performed cell communication analysis in all samples and each sample separately to explore robust signal pathways in the GBM cancer microenvironment. Through subgroup analysis for each included sample, we discovered that genes from the GALECTIN (LGALS9, CD44, CD45, and HAVCR2), GRN (GRN and SORT1), MK (MDK, PTPRZ1, NCL, ITGA6, ITGB1, LRP1, and SDC4), PTN (PTN, PTPRZ1, NCL, SDC3, and SDC4), SPP1 (SPP1, CD44, ITGAV, and ITGB1), ANNEXIN (ANXA1 and FPR1), and VISFATIN (NAMPT and INSR) signaling pathway networks were detected in most samples with consistent network patterns (; ; ). Specifically, M1 and M2 macrophages targeted other cells in the GALECTIN, GRN and SPP1 signaling pathway networks (), and GBM cancer cells regulated other cells in the MK and PTN signaling pathway network (; ).

Next, we conducted a Pearson correlation analysis to investigate whether m6A regulators were involved in the regulation of these signaling pathway networks (-). The GALECTIN signaling network plays an important role increasing induced regulatory T (iTreg) cell stability and suppressing T-cell proliferation (21). The expression level of genes in the GALECTIN signaling pathway network were higher in M1 macrophages, M2 macrophages, and T/NK cells (). Our Pearson correlation analysis showed that in M2 macrophages, expression levels of YTHDC1, YTHDC2, IGF2BP2, IGF2BP2, HNRNPC, ELAVL1, LRPPRC, ZC3H13, and ALKBH5 were positively correlated with expression of LGALS9, the gene encoding galectin-9 (). In T/NK cells, CD44 expression was positively correlated with that of YTHDF2 and HNRNPA2B1 and negatively correlated with FTO and ALKBH5. In T/NK cells, expression of HAVCR2, which encodes TIM-3, was positively correlated with that of YTHDF2. These results suggested that m6A modification may promote immunosuppressive activities through the GALECTIN signaling pathway network in the GBM immune microenvironment.

MDK Induces Macrophage Migration and M2 Polarization

Research reported that LRP1 attenuates proinflammatory macrophage activation as receptor for MDK (22). And our prementioned cell communication analysis found MK signaling pathway network (MDK, PTPRZ1, NCL, ITGA6, ITGB1, LRP1, and SDC4) in GBM samples. We assumed that the MK network (MDK/LRP1) played a significant role regulating macrophage activation in GBM. First, we collected four samples (Normal = 2; GBM = 2). IHC analysis showed a trend that the expression of MDK and LRP1 was higher in GBM patients (). To further validate our hypothesis, we first knockdown the expression of MDK in U87MG cell line with siRNA (), and coculture its supernate with induced THP-1. Tanswell assay revealed that compared with negative control (siCon), the knockdown of MDK (siMDK) significantly deceased the migration of macrophages (). Then, we induced THP-1 (human monocyte line) into macrophages in the upper chamber of transwell unit, and added recombinant protein MDK into the lower chamber of it. Interestingly, the addition of MDK protein showed significant effect on macrophages migration as well (). To assess the effect of MDK on the polarization of macrophages, we cocultured induced THP-1 with supernate from GBM in siMDK and siCon groups respectively. Flow cytometry showed less M2 polarization (CD11b+/CD206+) in the siMDK group, compared with siCon group (). The results suggested that MDK significantly induces an immunosuppressive macrophage differentiation. Finally, we conducted immunofluorescence staining to explore the co-localization of MDK and CD206 (marker of M2 macrophages). Results revealed more co-localization of MDK and CD206 in GBM tissues compared with that in normal brain tissues ().

Figure 5

MDK promotes migration and immunosuppressive polarization of macrophage. (A) The expression of MDK and LRP1 in GBM tissues (n = 2) and normal tissues (n = 2) determined using IHC (scale bar = 50 μm). And the CT and MRI images of GBM patients. (B) The siRNA knockdown effect confirmed with western blot and qPCR experiments. (C) THP-1 was treated with 185 ng/ml PMA for 24 h in the upper chamber of transwell unit. And then, coculture it with supernate from U87MG medium in the lower chamber for 24 h. The image showed the macrophages migrate through the membrane of chamber (scale bar = 100 μm). (D) THP-1 was treated with 185 ng/ml PMA for 24 h in the upper chamber of transwell unit. And then, coculture it with MDK protein at various concentration gradients in the lower chamber for 24 h. The image showed the macrophages migrate through the membrane of chamber (scale bar = 100 μm). (E) THP-1 was treated with 185 ng/ml PMA for 24 h. And then, coculture it with supernate from U87MG medium for 72 h The flow cytometry showed the expression levels of CD11b and CD206 in these cells. Bar plot showed the proportion of M2 macrophages (CD11b+/CD206+). **p < 0.01, ***p < 0.001, ****p < 0.0001. (F) Immunofluorescence staining showed the expression levels and co-localization of MDK and CD206 in normal brain tissues and GBM tissues.

m6A Regulators Correlate With ICPs in Immune Microenvironment

Having demonstrated that m6A regulators were associated with immune related process, we next addressed whether m6A regulators were correlated with ICPs in the GBM microenvironment (). Our results revealed that in GBM cancer cells, expression of CD274 (PD-L1), PDCD1LG2 (PD-L2), LGALS9 (Galectin-9), and PVR (CD155) were positively correlated with the expression of multiple m6A regulators, especially YTHDF2, YTHDF3, LRPPRC, METTL3, RBM15B, FTO, and ALKBH5 (). In M2 macrophages, LGALS9, CD86, and PVR showed positive correlations with YTHDC1, YTHDC2, YTHDF2, YTHDF3, ELAVL1, LRPPRC, METTL3, METTL14, ZC3H13, and ALKBH5 (). Also, we found positive correlations between HAVCR2 (TIM-3) and YTHDF2 () in T/NK cells. These results indicated that m6A modification may upregulate these suppressive immune check points and further promote the immunosuppressive microenvironment.

Bulk RNA-seq Analysis Identifies Two m6A-ICP Expression Patterns

As we have determined that ICPs in the GBM microenvironment were correlated with m6A regulators, we next explored the m6A-ICP expression pattern in bulk RNA-seq profiles. According to the similarities of expression of 23 m6A regulators and 13 ICPs, we performed consensus clustering, an unsupervised clustering algorithm, in TCGA, CGGA and REMBRANDT datasets, respectively. The optimal clustering stability was obtained when K = 2 (; ), and samples were clustered into two subgroups with different distinct features. Cluster 1 showed a low m6A/high ICP expression pattern in the heatmap of expression levels. A high m6A/low ICP expression pattern was observed for cluster 2 (; ). PCA for expression levels of m6A regulators and ICPs revealed prominent differences between these 2 subgroups (; ). Notably, the predicted potential therapeutic response of ICP inhibitors, conducted with the TIDE algorithm, was significantly different between these 2 subgroups in TCGA and REMBRANDT datasets (; P < 0.001). Patients in cluster 1 were more likely to response to ICP inhibitor therapy (> 50% were predicted responders), but fewer than 25% patients in cluster 2 were predicted to response to this immunotherapy. In the CGGA datasets, although no statistical significance was observed between these two clusters, cluster 1 showed a trend with higher response proportion to ICP inhibitor therapy (). However, the overall survival and clinical characteristics, including age, gender, IDH mutation status, and X1p19q codeletion status, were similar between these two clusters ().We further assessed the potential difference of functional states between these 2 clusters. An ssGSEA analysis detected consistent trends between TCGA and REMBRAMDT datasets with statistical significance in the following functional states: apoptosis, cell cycle, DNA damage, DNA repair, inflammation, quiescence and m6A modification (; ). In order to explore the immune infiltration of these 2 clusters, we utilized the CIBERSORT algorithm and compared the enrichment scores of 22 types of immune cells between these 2 clusters. Interestingly, in both TCGA and REMBRAMDT datasets, the proportion of follicular helper T cells in cluster 1 was lower than that in cluster 2, and there were more M2 macrophages in cluster 1 compared with cluster 2 (; ).

Figure 6

Bulk RNA-seq analysis for GBM patients and predictive model construction and validation. (A, B) Heatmap of m6A regulators and immune checkpoints (ICPs) revealed the different m6A-ICP expression patterns for 2 clusters identified by consensus clustering (for TCGA (A) and REMBRANDT (B) datasets separately). (C, E) Principal component analysis of 2 identified clusters (for TCGA (C) and REMBRANDT (E) datasets separately). (D, F) Predicted potential therapeutic response of ICP inhibitors of 2 identified clusters (for TCGA (D) and REMBRANDT (F) datasets separately). (G) Bar charts illustrating the differences of functional state scores between 2 identified clusters for TCGA dataset. (H) Bar charts illustrating the differences of CIBERSORT scores between 2 identified clusters for TCGA dataset. ns: P > 0.05, *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, ****P ≤ 0.0001. (I, J) Correlation network of the expression levels and methylation levels of m6A regulators and ICPs in 2 identified clusters. (K) The nomogram for distinguishing 2 identified clusters. (L, M) ROC curves of the nomogram distinguishing 2 identified clusters (L for training dataset and G for validation dataset). (N, O) Calibrate plots of the nomogram distinguishing 2 identified clusters (N for training dataset and O for validation dataset).

In revealing differences between the two clusters at the gene level, we detected 224 up-regulated DEGs and 33 down-regulated DEGs in cluster 1 compared with cluster 2 in TCGA dataset (). However, no DEGs were found in the REMBRAMDT dataset. Enrichment analyses of these 224 up-regulated DEGs uncovered immune-related BP terms (). KEGG analysis on the up-regulated DEGs was also conducted (). We further performed PPI network analyses for all DEGs, and the “cytoHubba” plugin from Cytoscape software was used to reveal hub genes. CASR, CCL19, CCR8, CCL13, RLN3, TAS2R42, SSTR4, MTNR1A, and RXFP3 were detected as hub genes ().

In order to address the potential regulatory mechanism between m6A regulators and ICPs at the multiomics level, samples with mRNA-seq profiles and DNA methylation profiles from TCGA were extracted for Pearson correlation analysis. Complex correlation networks of the expression levels and methylation levels of m6A regulators and ICPs for all samples and two clusters, respectively, were constructed (; ). Extensive correlation between expression levels of m6A regulators and ICPs were observed, and significant negative correlations between expression levels and DNA methylation levels were detected in IGF2BP3, METTLE3, ALKBH5, CTLA4, and TIGIT (). In the subnetworks of these 2 clusters, more positive correlations between expression levels of m6A regulators and ICPs were found in cluster 1. For IGF2BP3, a negative correlation between expression levels and DNA methylation levels were consistent in these two clusters. However, other statistically significant negative correlations between expression levels and DNA methylation levels were different between cluster 1 and 2 (). These results suggested different underlying mechanisms of the expression and DNA methylation of m6A regulators and ICPs between cluster 1 and 2 GBM patients.

Construction and Validation of Predictive Model

Having demonstrated that these 2 clusters had distinct features, we next constructed and validated a predictive nomogram model to distinguish these 2 clusters. The LASSO algorithm identified 32 genes, with the optimal λ being 0.02800602 (). Then, a best subset selection analysis was performed to determine the final model. Because of the limitations of the computing power of our device, we only analyzed models with up to 16 variables. The final model contains the following 6 variables: IGF2BP3, MAGEF1, PDCD1, CHKB, CD86, and TIGIT. The final multivariable logistic regression model was visualized as a nomogram (). The nomogram showed excellent discrimination performance in TCGA (training cohort, area under curve (AUC) = 0.991), REMBRANDT (validation cohort, AUC = 0.967) and CGGA (validation cohort, AUC = 0.721) datasets (; ). Calibration plots were drawn, and the calibration was as expected in both training and validation cohorts (; ). These evaluations proved that the nomogram model performed well in distinguishing these 2 GBM clusters.

Discussion

Treatment, especially immunotherapy, for GBM remains a challenge (1, 2). The role of m6A modification has been shown to be involved in cancer biology (4) and mediating immunotherapy resistance (6, 7). Therefore, there is an urgent need to apply new technologies, including single-cell analysis, to explore potential mechanisms of m6A modification in the GBM microenvironment. However, to date, no study has provided an m6A modification landscape for GBM microenvironment at the single-cell level. In this study, we retrieved single-cell RNA-seq datasets and used them to identify 7 types of cells in the GBM microenvironment and evaluated the associations between m6A regulators and functional states, potential BP, cell communication, and ICPs for these cell types. We discovered that m6A modification facilitates the stemness state in GBM cancer cells and promotes an immunosuppressive microenvironment through ICPs and the GALECTIN signaling pathway network. With bulk RNA-seq datasets, we further identified 2 clusters of patients with distinct m6A-ICP expression patterns (low m6A/high ICP and high m6A/low ICP) and predicted the responses to ICP inhibitors. A well performing nomogram model was constructed to distinguish these 2 GBM clusters.

The high inter-tumor heterogeneity of GBM immune microenvironment has reached a broad consensus. In the current research, we observed high inter-tumor heterogeneity of M1 and M2 macrophage proportions (). In around half of included samples, the infiltration proportions of M2 macrophage were higher than that of M1 macrophages. And higher proportions of M1 macrophages were found in the other samples (). Sørensen et al. (23) reported interesting phenomenon that on average 44% and 8% of the tumor associated macrophages (TAM) expressed the marker HLA-DR and TNF-α for M1 macrophages, respectively, while 10% and 3% of TAMs expressed IL10 and TGF-β1, which were markers for M2 macrophages. Also, these proportions varied greatly among different samples. These results indicated that the proportions of M1 and M2 macrophages in GBM have high inter-tumor heterogeneity.

There has been limited discussion about the distribution of expression of m6A regulators in different cell types and throughout cell cycle phases. The bulk RNA-seq dataset revealed that most m6A regulators have relatively high expression levels in GBM compared with levels in normal samples. This finding was especially apparent for HNRNPA2B1, HNRNPC and WTAP (). A previous study has shown increased expression of METTL3 and decreased expression levels of METTL14 and ALKBH5 and no significantly changes to the expression of FTO in GSCs (8). Another research detected elevated ALKBH5 in GSCs (24). The present study detected the expression levels of 23 m6A regulators in 7 types of cells and 3 cell cycle phases. We found that expression levels of m6A regulators were higher in M1 macrophages, GBM cancer cells, and T/NK cells than in other types of cells (), and cells at G1 phase showed higher expression levels of m6A regulators relative to that of cells at S and G2M phases (). HNRNPA2B1 and HNRNPC had obviously higher expression levels than other m6A regulators in all 7 types of cells and 3 cell cycle phases. Thus, besides the intratumor heterogeneity of m6A regulators, we also detected obvious intertumor heterogeneity of cell types and cell cycle phases in different samples (; ). These results suggested that the roles of m6A modification in GBM may be different in different types of cells, and further investigations of m6A modifications in GBM should take cell type into consideration.

Multiple studies have analyzed the mechanisms by which m6A modification maintains stemness of GSCs. For example, ALKBH5 promotes stemness of GSCs by sustaining FOXM1 expression (24, 25) and METTL3 ( 19), HNRNPA2B1 and HNRNPC ( 26) maintain stemness of GSCs by targeting SOX2. The product of the SOX2 gene is an oncogenic transcription factor in many cancers; accordingly, in the current research, we found significant positive correlations between SOX2 and YTHDC1, IGF2BP2, HNRNPC, ELAVL1, FMR1, KIAA1429, and ZC3H13 in GBM cancer cells (). Among these m6A regulators, IGF2BP2 had been proven to be a downstream m6A reader of SOX2 in colorectal carcinoma (20). Our results suggested that YTHDC1, IGF2BP2, HNRNPC, ELAVL1, FMR1, KIAA1429, and ZC3H13 may have potential interactions with SOX2 and further promote stemness of GBM cancer cells. The research of Su et al. demonstrated that FTO plays a carcinogenic role maintaining the self-renewal ability of GSCs via FTO/m6A/MYC/CEBPA signaling (27). Importantly, MYC is a commonly activated oncogene in human cancer. Our correlation analysis uncovered positive correlations between MYC and IGF2BP1, YTHDF1, YTHDF2, YTHDF3, HNRNPC, LRPPRC, METTL3, WTAP, KIAA1429, RBM15, RBM15B, CBLL1, FTO, and ALKBH5 (). The IGF2BP1 protein has been proposed to protect MYC mRNA from endonucleolytic attack (28). These results indicated a more extensive potential stemness regulation network between MYC and these m6A regulators in GBM cancer cells.

The m6A modification was reported to play a role in programmed cell death, including apoptosis (9, 29–32), ferroptosis (33, 34) and pyroptosis (35, 36). These studies mainly focused on cancer cells, and none of them have discussed the immune cells in the tumor microenvironment. In our research, we found that the m6A score was significantly correlated with apoptosis, ferroptosis and pyroptosis in M1 macrophages and apoptosis and pyroptosis in T/NK cells in GBM microenvironment (). Further correlation analysis of m6A regulators and these programmed cell death-related gene sets were visualized as regulation networks (). These results indicated a potentially complex mechanism of m6A modification-mediated programmed cell death in M1 macrophages and T/NK cells. It remains unknown whether programmed cell death in these immune cells could lead to the immunosuppressive microenvironment in GBM or not. And the regulatory role of m6A regulators in these processes is still unclear. Further research is needed in order to verify these findings and to analyze its impact on the GBM microenvironment.

The GALECTIN signaling pathway network includes 3 ligand-receptor pairs (LGALS9-CD44, LGALS9-CD45, LGALS9-HAVCR2; ). LGALS9 (Galectin-9) and HAVCR2 (TIM-3) are known as immune check points. CD45 (PTPRC) was an essential regulator of T cell activation (37), and CD44 is involved in diverse cellular processes including proliferation, apoptosis, and angiogenesis (38). Limited studies have discussed the role of m6A modification in the GALECTIN network. Wang et al. reported that IGF2BP1 stabilizes mRNA transcripts of CD44 and further promotes proliferation and invasion of GBM (39). A study from Lin et al. found positive correlations between expression of YTHDF2 and HAVCR2 in lower-grade glioma (40). In the current research, we observed robust occurrence of the GALECTIN network in included samples (), which indicated that the GALECTIN network plays an important role in the GBM microenvironment. The network mainly started from M1 and M2 macrophages and was targeted at other types of cells in the GBM microenvironment. Relatively high expression levels of genes in the GALECTIN network were detected in M1 macrophages, M2 macrophages, and T/NK cells (). Further correlation analysis found extensive significant correlations between m6A regulators and genes in the GALECTIN network (; ). These results indicated that the m6A modification was involved in the potential regulatory mechanism of the GALECTIN network. However, the detailed regulatory pathways between m6A modification and the GALECTIN network is still known. And whether inhibiting the GALECTIN network by regulation of m6A modification could lead to better prognosis of GBM patients remains to be further verified.

We conducted a series of in vitro experiments to validate the regulation role of MDK in GBM microenvironment (). Based on our results, MDK, which was secreted by GBM cancer cells, induced the migration and immunosuppressive polarization of macrophages. Our results were consistent with the role of MDK identified by Zhang et al. in 2021 (41). In their gallbladder cancer research, MDK interacts with its receptor LRP1, which is expressed by tumor-infiltrating macrophages, and further promotes immunosuppressive macrophage differentiation. We assume that MDK could be a treatment target for GBM, further in vivo validation and translational research are expected to develop a novel therapy for GBM.

Rapid progress in cancer immunotherapy, especially in the development of ICP inhibitors, has revolutionized the treatment of many solid tumors and driven the study of immunotherapy in glioma. However, the effect of ICP blockades in GBM has not been satisfactory (42). Some studies have found that m6A modification was correlated with immune infiltration in glioma (43–47). However, to date, no research has analyzed the relationship between m6A modification and ICPs in GBM. In the current study, we revealed extensive correlations between m6A regulators and ICPs in GBM cancer cells, M2 macrophages, and T/NK cells (). Notably, in GBM cancer cells, CD274 (PD-L1), PDCD1LG2 (PD-L2), LGALS9 (Galectin-9), and PVR (CD155) were correlated with m6A regulators. In the bulk RNA-seq analysis, we identified 2 clusters of patients. Cluster 1 patients (low m6A/high ICP) showed higher predicted response rates to ICP inhibitors, and cluster 2 patients (high m6A/low ICP) had lower predicted response rates (). Although these 2 clusters showed different m6A-ICP expression patterns, they had similar overall survival times (). These results indicated that these 2 clusters may have different mechanisms for immune escape and require different treatment strategies.

In the current research, the LASSO analysis and best subset selection analysis filtered 6 genes, which could distinguish these 2 clusters. Among them, IGF2BP3 was a m6A reader. PDCD1, CD86, and TIGIT were significant ICPs. MAGEF1 (melanoma antigen family F1) is a member of MAGE family and belongs to type 2 MAGE (T2M) category (48). Arora et al. ’s research found that in glioma downregulation of T2Ms was associated with immune infiltration and poor overall survival (48). Weon et al. reported that MAGEF1 alters DNA repair enzymes via the cytosolic iron-sulfur cluster assembly (CIA) pathway (49). And MAGEF1 was highly amplified in multiple human cancer types and related with increased mutational burden (49). Based on these findings, MAGEF1 might related with ICPs via m6A modification or CIA pathways, which needed to be further explored. CHKB (choline kinase beta) played a key role maintaining the normal phosphatidylcholine level (50). The relationship between CHKB and glioma haven’t been explored. Our results showed that lower expression level of CHKB suggested higher possibility of cluster 2 (high m6A/low ICP; ). However, the potential mechanisms need to be further analyzed.

Based on the correlations between m6A regulators and ICPs found in the current study, we suggest that by targeting these related m6A regulators, researchers may downregulate the expression of these immune suppressive ICPs and further enhance the efficacies of ICP inhibitors for cluster 2 patients. Further research is needed to verify these findings and to specify how it influences strategies directing the use of ICP inhibitors for GBM.

Research have found that tumors with extensive infiltration of immunosuppressive macrophages are refractory to ICP inhibitor therapy (51, 52). However, contrary results were observed in our research. TIDE algorithm predicted that cluster 1 was more responsive to ICP inhibitors (), and CIBERSORT algorithm showed higher M2 macrophage infiltration in cluster 1(; ). These results may reflect the inner complexity of GBM. Further high-quality research of ICP inhibitors in GBM are expected.

The current research revealed the m6A regulators expression landscape of GBM at single-cell level. Our research suggested that by altering the m6A modification in GBM, researchers may be able to influence the stemness status and immunosuppressive microenvironment of GBM. Combined with immunotherapy, these regulations have potential to advance the treatment effect of immunotherapy.

Some limitations in the current study should be acknowledged. First, as a retrospective study, these analyses were conducted using published datasets, and the results about potential mechanisms of action of m6A modifications in the GBM microenvironment need further experimental validation. Second, we were not able to identify clusters of GSCs and to separate the T/NK cluster into more detailed clusters. Here, the problem may be due to limited cell numbers for these cell types. Third, missing data was observed after merging the 5 single-cell datasets. However, the data of m6A regulators was almost complete, and only a few datasets partially lacked genes for some ICPs (). Also, the information regarding the rest of the genes is complete. To avoid the influence of the missing data in the analysis, we excluded samples with missing data in the analyses of correlations between m6A regulators and ICPs ().

Conclusion

Through analyses at the single-cell level, for the first time, we discovered that m6A modification facilitates stemness state in GBM cancer cells and promotes the immunosuppressive microenvironment through ICPs and the GALECTIN signaling pathway network in the GBM microenvironment. We further identified 2 clusters of patients with distinct m6A-ICP expression patterns (low m6A/high ICP and high m6A/low ICP) and predicted the response of ICP inhibitors in bulk RNA-seq analysis. A well performing nomogram model was constructed to distinguish these 2 GBM clusters. We hope the novel understanding of the roles of m6A modification in the GBM microenvironment may assist the development of immunotherapy and precision treatment in GBM.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: The data analyzed in this study were obtained from Gene Expression Omnibus (GEO) at GSE141383, GSE138794, GSE84465, GSE103224, GSE89567, GSE124731 and GSE108474 (REMBRANDT), and from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/).

Ethics Statement

The studies involving human participants were reviewed and approved by Jinling Hospital research ethics board. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

CM, XC and SZ were involved in conception and design of the study. FY, XC, YG, YA, CD, JZ, JY, CT, and ZC were involved in development of methodology, acquisition of data, analysis and interpretation of data. XC was involved in writing of the manuscript. YW and AZ contributed to the interpretation of the experimental results. CM oversaw the conception and design of the study, interpretation of data and writing of the manuscript. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.