Long Noncoding RNA IRAIN Drives Immunometabolic Reprogramming in Glioma via the IGF1R-JAK-STAT-BIRC5 Axis

Glioma represents the most common primary malignancy of the central nervous system and remains among the deadliest human cancers despite progress in neurosurgery, radiotherapy, and chemotherapy. High-grade gliomas, particularly glioblastoma multiforme (GBM), are characterized by rapid proliferation, diffuse infiltration, and inevitable recurrence^1,2. Recently, a range of effective therapeutic approaches has arisen, such as surgical interventions, immunotherapy, and chemotherpy³. Median survival remains approximately 14-18 months, even with maximal therapy, underscoring the need for novel therapeutic strategies and prognostic biomarkers⁴.

The interplay between the tumor microenvironment (TME) and cancer cells critically drives the aggressive growth and molecular heterogeneity of gliomas. Numerous immune constituents, including microglia, neuronal precursor cells, vascular cells, and adaptive immune cells, are intricately involved in shaping the TME^5,6,7. Recent years have witnessed remarkable advances in immunotherapy, ushering in a transformative era in oncology. Nevertheless, immune checkpoint blockade (ICB) and other immunotherapeutic approaches have shown only limited success in glioma⁸. The failure of immunotherapy in this context stems from several intrinsic features of the glioma microenvironment-namely, low tumor mutational burden, pronounced intratumoral heterogeneity, poor blood-brain barrier permeability, downregulation of major histocompatibility complex molecules, and sparse T-cell infiltration^9,10,11,12.

Increasing evidence suggests that the tumor microenvironment (TME), which exhibits metabolic dysfunction, is marked by acidic conditions, nutrient deprivation, and the accumulation of metabolites that suppress immune responses. This condition subsequently leads to the dysfunction of immune cells that infiltrate tumors and diminishes the effectiveness of immunotherapy¹³. Simultaneously, immune cells within the TME undergo context-dependent metabolic shifts that dictate their differentiation and effector functions. For instance, activated T cells rely on glycolysis for rapid proliferation, whereas regulatory T cells (Tregs) and tumor-associated macrophages (TAMs), particularly the M2 subtype, depend on oxidative phosphorylation and fatty-acid metabolism to maintain their immunosuppressive activity^14,15.

Therefore, integrating immune and metabolic signatures offers a rational framework to better classify gliomas and predict patient prognosis. Long noncoding RNAs (lncRNAs) have recently emerged as versatile regulators of chromatin architecture, gene transcription, and post-transcriptional modification. Several lncRNAs have been implicated in glioma pathogenesis, influencing proliferation, invasion, and immune modulation¹⁶. Among them, the imprinted antisense transcript IRAIN (IGF1R antisense imprinted non-protein coding RNA) has garnered attention for its role in regulating the insulin-like growth factor 1 receptor (IGF1R) gene through cis-acting chromatin interactions, which has been proven to be downregulated and decreased cis competition control, especially in cancer cells such as breast cancer, leukemia cell lines, laryngeal cancer, and acute myeloid leukemia, which directly leads to upregulation of IGF1R^17,18,19. Because IGF1R signaling activates downstream JAK-STAT and PI3K-AKT pathways -- both central to tumor growth and immune regulation^20,21 -- alterations in IRAIN expression may critically affect tumor immunometabolism. However, the role of IRAIN in glioma and its potential impact on the immunometabolic landscape have not yet been elucidated. Addressing this knowledge gap could provide mechanistic insight into how noncoding RNAs coordinate metabolic adaptation and immune evasion in glioma.

In this study, we combined multi-cohort transcriptomic data with experimental validation to explore immunometabolic regulation in glioma. Using integrated machine-learning frameworks across the TCGA, CGGA (693/325), and GEO (GSE43378) cohorts, we identified a robust 17-gene immunometabolic-related gene (IMRG) prognostic signature capable of stratifying patient survival across datasets. We subsequently focused on IRAIN as a potential upstream regulator of this network, given its predicted interaction with the IGF1R-JAK-STAT signaling axis. Functional assays revealed that IRAIN overexpression suppresses glioma proliferation and angiogenesis by inhibiting IGF1R expression and downstream STAT3 phosphorylation, resulting in reduced expression of the anti-apoptotic protein BIRC5 (also known as Survivin). Taken together, these findings support a working model in which IRAIN functions as a mechanistic bridge linking immune regulation and metabolic reprogramming in glioma. We hypothesize that dysregulated IRAIN expression promotes tumor progression by disturbing this balance, thereby fostering an immunosuppressive microenvironment and conferring resistance to therapy.

All procedures involving human tissues complied with institutional guidelines and the Declaration of Helsinki and were approved by the Institutional Review Board of Fujian Medical University (Approval No. 2021KYB089). Written informed consent was obtained from all participants prior to tissue procurement.

Gene expression and survival analysis
RNA sequencing data and corresponding clinical information were obtained from multiple public databases. 1) TCGA cohort: RNA-seq (FPKM) data for 175 glioblastoma multiforme (GBM) and 534 low-grade glioma (LGG) samples were downloaded from The Cancer Genome Atlas (https://portal.gdc.cancer.gov/); 2) Normal controls: Expression profiles of 211 normal brain tissues and 662 glioma tissues were downloaded from the UCSC Xena database (https://xenabrowser.net/datapages/); 3) External validation: Data from CGGA693 and CGGA325 cohorts were obtained from the Chinese Glioma Genome Atlas (http://www.cgga.org.cn); 4) GEO dataset: The GSE43378 dataset, containing expression and clinical data for 50 glioma samples, was downloaded from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). All raw count data were converted to transcripts per million (TPM) and log₂-transformed. For datasets already normalized, expression matrices were examined to ensure comparable distributions. Genes with TPM values < 1 in more than 80% of samples were excluded. Missing clinical information (age, IDH status, 1p/19q codeletion, MGMT methylation) was removed using complete-case filtering. Batch effects among datasets were adjusted using the ComBat algorithm implemented in the R package sva. Expression values were standardized by z-score transformation within each dataset. Survival analyses were performed using the R packages survival and survminer. Patients were dichotomized into high- and low-expression groups according to the median expression level of IRAIN. Kaplan-Meier survival curves were generated, and statistical significance was evaluated by the log-rank test. Hazard ratios (HRs) and 95% confidence intervals (CIs) were estimated using Cox proportional hazards regression models.

Definition of immune and metabolic gene sets
Immune-related genes (IRGs, n = 2,483) were obtained from the ImmPort database (https://www.immport.org/shared/), and metabolic-related genes (MRGs, n = 948) were obtained from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/). The combined set of these genes was defined as immunometabolic-related genes (IMRGs). These gene lists served as references for subsequent differential expression and network analyses.

Differential expression and weighted gene co-expression network analysis
Differentially expressed genes (DEGs) between normal brain and glioma tissues were identified using the R package limma. Expression data were fitted with a linear model followed by empirical Bayes moderation. Genes with |log₂ fold change| > 1.5 and false discovery rate (FDR) < 0.05 were considered significantly differentially expressed. Weighted gene co-expression network analysis (WGCNA) was conducted using the R package WGCNA. Outlier samples were excluded through hierarchical clustering. The soft-thresholding power was set to β = 8 to achieve a scale-free topology fit index (R² ≥ 0.85) while maintaining adequate mean connectivity. Topological overlap matrices (TOM) were constructed, and genes were grouped into modules with a minimum size of 50 using the dynamic tree-cut algorithm. Module eigengenes were correlated with clinical traits, and the module most strongly associated with glioma (Pearson's r > 0.7, P < 1×10^-10) was selected for hub gene identification.

Machine learning-based prognostic model construction
A comprehensive leave-one-out cross-validation (LOOCV) framework integrating ten machine-learning algorithms was applied to construct and evaluate prognostic models. In total, 101 combinatorial workflows were implemented using the TCGA cohort as the training dataset. Prognosis-associated immunometabolic-related genes (IMRGs) were first identified by univariate Cox regression (P < 0.05). The optimal model was determined by maximizing the mean Harrell's concordance index (C-index) across three validation datasets (CGGA693, CGGA325, and GSE43378). The resulting RSF-Enet model (α = 0.3) demonstrated the highest predictive performance and maintained robust generalizability across independent cohorts²².

TME and immune infiltration
To comprehensively characterize the immunogenomic landscape, we employed a multi-tiered analytical approach. First, immune and stromal infiltration levels were quantified using the ESTIMATE algorithm²³. Differential expression of key immune checkpoint molecules, including PDCD1, CTLA4, and LAG3, was then assessed through limma-based analysis, and the correlations among checkpoint genes were visualized using correlation matrices. Somatic mutation profiles from 903 glioma samples in the TCGA cohort were used to calculate tumor mutational burden (TMB), microsatellite instability (MSI), and tumor immune dysfunction and exclusion (TIDE) scores to predict potential responses to immunotherapy. Patients were subsequently stratified into four prognostic groups according to combined TMB status (high/low) and risk scores (high/low), and survival outcomes were compared using Kaplan-Meier analysis.

Functional enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted using the R package clusterProfiler. Enrichment results with adjusted P values < 0.05 were considered statistically significant. Overrepresented biological processes, cellular components, and molecular functions were visualized using dot plots and bar plots. Protein-protein interaction (PPI) networks were constructed using the STRING database (≥ 0.4) and visualized in Cytoscape. Functional modules within the PPI network were identified using the MCODE algorithm. Gene-gene interaction and co-expression networks were further analyzed using GeneMANIA (https://string-db.org; confidence score ≥ 0.4) and visualized in Cytoscape. Functional modules within the PPI network were identified using the MCODE algorithm. Gene–gene interaction and co-expression networks were further analyzed using GeneMANIA (https://genemania.org), which integrates information on physical and genetic interactions, shared pathways, and co-expression patterns to infer potential functional associations.

Clinical specimens
Fresh glioma tissues (n = 6) and paired adjacent non-tumorous brain tissues (n = 6; located at least 3 cm from the tumor margin and histologically confirmed as tumor-free) were collected from patients undergoing primary glioma resection at the Zhangzhou Affiliated Hospital of Fujian Medical University. None of the patients had received chemotherapy or radiotherapy prior to surgery. All pathological diagnoses were independently verified by two neuropathologists according to the 2021 World Health Organization (WHO) classification of central nervous system tumors. Immediately after surgical excision, tissue specimens were rinsed with ice-cold phosphate-buffered saline (PBS) to remove residual blood, snap-frozen in liquid nitrogen (-196 °C), and stored at -80 °C until RNA extraction.

Cell lines and cell culture
Human glioblastoma cell lines SHG44, U251, A172, and T98G, as well as normal human glial cells (HEB), were obtained from authenticated repositories and confirmed to be free of mycoplasma contamination prior to use. Cells were maintained in Dulbecco's Modified Eagle's Medium (DMEM, high glucose) supplemented with 10% fetal bovine serum (FBS), 2 mM L-glutamine, and 1% penicillin-streptomycin, at 37 °C in a humidified incubator with 5% CO₂. Cells were passaged every 4-5 days upon reaching 80-90% confluence. To establish IRAIN-overexpressing and control cell lines, cells were transduced with lentiviral vectors carrying the full-length IRAIN transcript or an empty vector as control. Stable clones were selected using puromycin (2 µg/mL) for 14 days. Overexpression efficiency was confirmed by quantitative reverse transcription PCR (qRT-PCR) prior to downstream assays.

3- (4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) cell proliferation assay
Cells were seeded in 96-well plates at a density of 1 × 10⁴ cells per well in 100 µL of complete culture medium. At 24, 48, and 72 h after seeding, 20 µL of MTT solution (5 mg/mL in phosphate-buffered saline) was added to each well and incubated for 4 h at 37 °C. The supernatant was then removed, and 150 µL of dimethyl sulfoxide (DMSO) was added to dissolve the formazan crystals. The plate was gently agitated for 10 min to ensure complete solubilization. Absorbance was measured at 490 nm using a microplate spectrophotometer. Background readings from blank wells were subtracted. Cell viability was calculated relative to the 24-h or control group (set as 1.0). All experiments were performed with six technical replicates and three independent biological replicates. Data are expressed as mean ± standard deviation (SD), and statistical significance was determined using a two-tailed t-test.

Flow cytometry for apoptosis (Annexin V - FITC/PI staining)
Cells were seeded at 60-70% confluence and treated for 24 h under the indicated conditions. Floating and adherent cells were collected using EDTA-free trypsin, combined, and washed twice with ice-cold PBS. Cell pellets were resuspended in Annexin V binding buffer (10 mM HEPES pH 7.4, 140 mM NaCl, 2.5 mM CaCl₂) at 1 × 10⁶ cells/mL. For each sample, 100 µL of suspension was incubated with 5 µL of Annexin V-FITC and 5 µL of propidium iodide (PI; 50 µg/mL stock) in the dark for 15 min at room temperature. Following the addition of 400 µL of binding buffer, samples were kept on ice and analyzed within 1 h on a flow cytometer (488 nm excitation; 530/30 nm for FITC and >585 nm for PI). Appropriate single-stain and fluorescence-minus-one controls were included for compensation. At least 10,000 events were recorded per sample. Data were analyzed by quadrant gating: live (Annexin V⁻/PI⁻), early apoptotic (Annexin V⁺/PI⁻), late apoptotic (Annexin V⁺/PI⁺), and necrotic (Annexin V⁻/PI⁺) populations. Percentages of early + late apoptotic cells were reported (mean ± SD, n = 3).

Quantitative real-time PCR (qRT-PCR)
Total RNA was isolated using an acid phenol-guanidinium reagent according to the manufacturer's protocol. RNA purity was verified by spectrophotometry (A₂₆₀/A₂₈₀ = 1.8-2.1), and integrity was confirmed by gel electrophoresis (RNA integrity number ≥ 7). One microgram of total RNA was treated with DNase I and reverse-transcribed in a 20 µL reaction using random hexamers and oligo(dT) primers. The reaction was carried out at 25 °C for 10 min, 50 °C for 30 min, and 85 °C for 5 min. Quantitative PCR was performed in a 10 µL system containing 5 µL of 2× SYBR Green Master Mix, 0.3 µM each primer, and 1 µL of cDNA (≈ 20 ng RNA equivalent). Thermal cycling conditions were 95 °C for 5 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 30 s, then a melt-curve analysis from 65 °C to 95 °C in 0.3 °C increments. All reactions were performed in triplicate, together with no-template and minus-RT controls. Ct values > 35 or technical replicate SD > 0.5 were excluded. Relative expression was calculated using the 2⁻ΔΔCt method, with GAPDH as the internal control. Mean ± SD values from three independent biological replicates were reported, and group differences were analyzed using a two-tailed t-test.

Western blot analysis
Cells were lysed on ice in RIPA buffer (50 mM Tris-HCl, pH 7.4, 150 mM NaCl, 1% NP-40, 0.5% sodium deoxycholate, 0.1% SDS) supplemented with protease and phosphatase inhibitors. Lysates were incubated for 30 min on ice with intermittent vortexing and cleared by centrifugation at 12,000 × g for 15 min at 4 °C. Protein concentrations were measured by BCA assay, adjusted to 1-2 µg/µL, and mixed 1:3 with 4× Laemmli buffer (final 1× buffer containing 100 mM DTT). Samples were denatured at 95 °C for 5 min. Equal amounts of protein (50 µg) were resolved by 12 % SDS-PAGE at 100 V for 90 min and electro-transferred to PVDF membranes at 250 mA for 90 min. Membranes were blocked with 5 % non-fat milk in TBST (0.1 % Tween-20) for 1 h at room temperature (or 5 % BSA for phosphoproteins) and incubated overnight at 4 °C with primary antibodies against IGF1, IGF1R, JAK2, p-JAK2 (Y1007/1008), STAT3, p-STAT3 (Y705), BIRC5, and β-actin (typical dilution 1:1000, β-actin 1:5000). After three 10-min washes in TBST, membranes were incubated with HRP-conjugated secondary antibody (1:5000) for 1 h at room temperature, washed again, and developed using chemiluminescent substrate. Band intensities were quantified with ImageJ, normalized to β-actin or total protein, and expressed as mean ± SD from three independent experiments.

Immunocytochemistry
Cells grown on sterile glass coverslips were rinsed twice with PBS and fixed in 4% paraformaldehyde for 15 min at room temperature. After three PBS washes, cells were permeabilized with 0.2% Triton X-100 for 10 min, blocked with 5% bovine serum albumin (BSA) for 1 h, and incubated overnight at 4 °C with primary anti-CD31 antibody (1:200 dilution in 1% BSA). Following three PBS washes, cells were incubated with Alexa Fluor-conjugated secondary antibody (1:500 dilution) for 1 h in the dark, counterstained with DAPI (1 µg/mL, 5 min), and mounted in antifade medium. Images were captured using a fluorescence microscope under identical exposure and gain settings. The percentage of CD31-positive area was quantified in five randomly selected non-overlapping fields per sample using ImageJ software. This assay was performed in cell models rather than tissue sections.

Statistical analysis
Statistical analyses were performed utilizing R version 4.3.0 along with its associated packages. To compare categorical variables, the chi-squared test was employed, while continuous variables were assessed using either the Wilcoxon rank-sum test or the T test. The evaluation of continuous variables was accomplished through Pearson's correlation coefficient. Survival analyses were carried out using the survival package, which included Cox proportional hazards modeling and the generation of Kaplan-Meier curves, with optimal stratification thresholds established by the survminer package and the formula Riskscore = . The CompareC package was used to assess the C-indices of various variables. The receiver operating characteristic curve (ROC), aimed at predicting binary categorical variables, was generated using the pROC package. Additionally, the time-dependent area under the ROC curve (AUC) for survival metrics was analyzed using the timeROC package. All statistical tests were conducted with a two-sided approach. A significance level of P < 0.05 was considered statistically significant.

IRAIN is downregulated in glioma and associated with adverse clinicopathologic features
The analysis process is illustrated in Figure 1. IRAIN expression was first assessed by qRT-PCR in primary astrocytes, glioma tissues and glioma cell lines. IRAIN levels were markedly reduced in both low-grade and high-grade glioma tissues and in all four glioma cell lines (SHG44, A172, U251 and T98G) compared with normal astrocytes (Figure 2A). Consistent with these findings, analysis of RNA-seq data from the TCGA and CGGA cohorts showed significantly lower IRAIN expression in glioma samples than in normal brain tissues (Figure 2B). When patients were stratified by IRAIN expression, those in the low-IRAIN group had significantly shorter overall survival and progression-free survival than those with high IRAIN expression (both P < 0.001; Figure 2C). Subgroup analyses further revealed that IRAIN expression was decreased in 1p/19q non-codeleted, IDH-wildtype, MGMT-unmethylated, and GBM cases, as well as in older patients (>65 years), whereas no significant difference was observed between male and female patients (Figure 2D). Gene Ontology and KEGG enrichment analyses of IRAIN-associated genes indicated significant enrichment in chromatin-related processes and cancer- and metabolism-related signaling pathways, including JAK-STAT, AMPK, mTOR, FoxO, and cellular senescence pathways (Figure 2E), suggesting that IRAIN loss is linked to malignant biological behavior and an altered immunometabolic program in glioma.

Identification of IRAIN-associated immunometabolic genes
To identify IRAIN-related genes with potential immunometabolic functions, we first performed differential expression analysis of TCGA glioma samples stratified by IRAIN expression. This analysis identified 4,420 differentially expressed genes (DEGs; |log2FC| > 1.5, FDR < 0.05; Figure 3A), including 3,292 upregulated and 1,128 downregulated genes in the low-IRAIN group. We next constructed a weighted gene co-expression network (WGCNA) to capture IRAIN-related co-expression patterns. A soft-thresholding power of β = 8 was selected because it yielded an approximate scale-free topology (signed R2 ≈ 0.85) while maintaining reasonable mean connectivity (Figure 3B). Among the resulting modules, the brown module showed the strongest positive correlation with the IRAIN-associated trait (r = 0.77, P = 2 × 10-170), whereas the cyan and midnight-blue modules were negatively correlated (Figure 3C, left). Within the brown module, module membership was strongly associated with gene significance for the trait (r = 0.73, P < 1 × 10-200; Figure 3C, right), and these genes were therefore considered IRAIN-related hub genes. Intersecting DEGs, brown-module hub genes, and curated immune/metabolism-related gene sets yielded 59 immunometabolic candidate genes (IMRGs; Figure 3D). Among these, univariate Cox regression analysis in the TCGA cohort identified 49 genes significantly associated with overall survival (P < 0.05), including 38 risk genes (HR > 1) and 11 protective genes (HR < 1; Figure 3E). These 49 prognostic IMRGs were used for subsequent model construction.

Construction and validation of a 17-gene immunometabolic risk signature
Using the 49 prognostic IMRGs, we applied a leave-one-out cross-validation framework that evaluated 101 combinations of 10 machine-learning algorithms across the TCGA training set and three external validation cohorts. As shown in Figure 4A, the Random Survival Forest (RSF) + Elastic Net model with α = 0.3 achieved the highest and most stable concordance indices across datasets and was therefore selected as the final model. This procedure yielded a 17-gene IMRG signature, with positive and negative coefficients indicating risk- and protection-associated genes, respectively (Figure 4B). Based on the 17-gene risk score, patients in the TCGA, CGGA693, CGGA325 and GSE43378 cohorts were stratified into high- and low-risk groups. In all four datasets, high-risk patients showed significantly shorter overall survival than low-risk patients (log-rank P < 0.001; Figure 4C-F). Time-dependent ROC curves confirmed good discriminative ability, with 1-, 3-, and 5-year AUCs of 0.893, 0.904, and 0.884 in the TCGA training cohort and comparable performance in the validation cohorts.

To determine whether the IMRG risk score provides prognostic information beyond conventional clinicopathologic factors, we integrated it with age, sex, WHO grade, IDH status, 1p/19q codeletion status, and MGMT promoter methylation. Univariate Cox regression identified grade, IDH status, 1p/19q status, MGMT methylation, and the IMRG risk score as significant predictors of overall survival (Figure 5A). In multivariate analysis, the risk score remained an independent prognostic factor after adjustment for these variables (P < 0.001; Figure 5B). Consistent with its association with aggressive disease, higher risk scores were observed in tumors with higher WHO grade, IDH wild-type, 1p/19q non-codeletion, MGMT unmethylation, and in older patients, whereas no obvious difference was seen between sexes (Figure 5C). Based on the IMRG risk score and clinical variables, we constructed a nomogram to estimate 1-, 3-, and 5-year overall survival (Figure 5D). The nomogram showed excellent performance, with a C-index of 0.879 (95% CI 0.855-0.903). Calibration plots indicated close agreement between predicted and observed survival, and ROC as well as decision-curve analyses demonstrated that the combined model outperformed individual clinical parameters. Across the TCGA, CGGA693, CGGA325, and GSE43378 cohorts, risk-score distributions, survival status plots, and heatmaps of the 17 IMRGs consistently showed that deaths were enriched in the high-risk groups and that these patients displayed distinct IMRG expression patterns (Figure 6A-D). Together, these findings validate the 17-gene IMRG signature as a stable, independent prognostic marker linking immunometabolic dysregulation to glioma progression.

Association of the IMRG risk signature with the tumor immune microenvironment
Comprehensive profiling of the tumor microenvironment (TME) revealed marked immune differences between IMRG-defined risk groups. Compared with low-risk tumors, high-risk gliomas exhibited significantly altered enrichment of multiple immune cell populations and immune-related functions, including antigen-presenting cells, T-cell subsets, and interferon responses (Figure 7A). When patients were stratified by IRAIN expression, numerous immune cell subsets also differed between IRAIN-high and IRAIN-low groups, indicating that IRAIN dysregulation is closely linked to immune infiltration patterns (Figure 7B). TME scores calculated by the ESTIMATE algorithm (stromal, immune, and ESTIMATE scores) varied significantly between high- and low-risk groups and between IRAIN-high and IRAIN-low tumors (Figure 7C), suggesting broad remodeling of both stromal and immune compartments. LinkET-based correlation analysis further revealed significant associations of IRAIN and representative IMRGs such as BIRC5 with specific immune cell subsets (Figure 7D). Consistent with an immunosuppressive milieu, most immune checkpoint genes -- including PDCD1, CTLA4, LAG3, HAVCR2, and SIGLEC15 -- were more highly expressed in the high-risk group (Figure 7E). In addition, high-risk patients had significantly higher TIDE scores than low-risk patients (Figure 7F), implying stronger immune evasion and potential resistance to immune checkpoint blockade. Together, these findings indicate that IRAIN downregulation and a high IMRG-based risk score are associated with a more immunosuppressive, therapy-resistant TME in glioma.

Somatic mutation landscape and functional enrichment of IMRG-related genes
Somatic mutation analysis was performed to explore the genomic background of the IMRG-based risk model. In the TCGA cohort, low- and high-risk groups displayed distinct mutational landscapes (Figure 8A, upper panels). Low-risk tumors were characterized by frequent IDH1, TP53, and ATRX mutations, a pattern typical of less aggressive gliomas, whereas high-risk tumors more commonly harbored PTEN, TTN, and EGFR alterations, which are associated with malignant progression. We further evaluated the prognostic impact of tumor mutational burden (TMB). Patients with low TMB had significantly longer overall survival than those with high TMB (P < 0.001; Figure 8A, lower left). When TMB status was combined with the IMRG risk score, the high-TMB/high-risk group showed the poorest outcome, while the low-TMB/low-risk group had the most favorable survival, with the other two groups in between (P < 0.001; Figure 8A, lower right). GeneMANIA analysis of the 17 IMRGs and closely related genes revealed a densely connected interaction network with EGFR at the core and extensive links to IDH1, BIRC5, CDKN2B, and several other nodes (Figure 8B). GO and KEGG enrichment analyses indicated that these genes were mainly involved in immune regulation and cell-cycle-related processes, including antigen processing and presentation, regulation of NF-κB signaling, glutathione metabolism, human cytomegalovirus infection, and the JAK-STAT signaling pathway (Figure 8C). STRING-based protein-protein interaction (PPI) analysis further highlighted central clusters organized around EGFR, BIRC5, RRM1, IDH1, BMPR1A, and CDKN2B (Figure 8D). These findings suggest that the IMRG risk signature is embedded in interconnected pathways related to cell-cycle control, redox metabolism, and inflammatory signaling, providing a mechanistic basis for subsequent experimental interrogation of IRAIN-JAK-STAT-BIRC5 crosstalk in glioma cells.

IRAIN suppresses glioma proliferation and angiogenesis via the IGF1/IGF1R-JAK2-STAT3-BIRC5 axis
We next examined the functional consequences of IRAIN overexpression in glioma cells. IRAIN exerted clear tumor-suppressive effects by inhibiting proliferation, promoting apoptosis and reducing angiogenic capacity. In CCK-8/MTT assays, IRAIN overexpression significantly decreased cell viability at 24-72 h compared with the control group; this growth-inhibitory effect was further enhanced by IGF1 knockdown and was partially rescued by SURVIVIN (BIRC5) overexpression (Figure 9A). Annexin V-FITC/PI flow cytometry confirmed that IRAIN markedly increased the proportion of apoptotic cells, with the highest apoptotic rate observed in the IRAIN + si-IGF1 group, whereas SURVIVIN overexpression attenuated IRAIN-induced apoptosis (Figure 9B). Immunocytochemical staining for CD31 showed a pronounced reduction of endothelial marker expression in IRAIN-overexpressing cells; this decrease was further strengthened by si-IGF1 and largely restored by SURVIVIN (Figure 9C), indicating impaired angiogenic potential. Consistently, western blotting demonstrated that IRAIN reduced the levels of IGF1, IGF1R, JAK2, STAT3, and BIRC5/SURVIVIN, whereas STAT3 overexpression reversed the downregulation of BIRC5 (Figure 9D). These results support a causal regulatory axis in which IRAIN attenuates IGF1/IGF1R signaling, thereby suppressing downstream JAK2-STAT3 activation and reducing expression of the anti-apoptotic effector BIRC5. Overall, IRAIN emerges as a central regulator linking immunometabolic signaling to glioma malignancy and represents a potential therapeutic target for restoring immune and metabolic homeostasis within the tumor microenvironment.

DATA AVAILABILITY:
The RNA sequencing (RNA-seq) and clinical data of 175 GBM and 534 LGG were downloaded from the Cancer Genome Atlas (TCGA) database. For external validation, gene expression data along with the corresponding clinical information of glioma patients were obtained from the Chinese Glioma Genome Atlas (CGGA) database and the GEO database. All data generated or analyzed during this study have been fully presented and discussed within the manuscript. The raw data of this study are uploaded as Supplementary File 1.

Figure 1: Overview of the analytical workflow for identifying IRAIN-associated immunometabolic regulatory genes in glioma. Schematic illustration of the study design integrating bioinformatics and experimental validation. Transcriptomic datasets (TCGA, CGGA, GEO) were analyzed to identify immunometabolic regulatory genes (IMRGs). Weighted gene co-expression network analysis (WGCNA), machine-learning model construction, and functional validation in glioma cell lines were performed to elucidate the mechanistic role of IRAIN in regulating the IGF1R-JAK-STAT-BIRC5 axis. Please click here to view a larger version of this figure.

Figure 2: IRAIN expression patterns and associated signaling pathways in glioma. (A) qRT-PCR analysis of IRAIN expression in primary astrocytes, low-grade and high-grade glioma tissues, and glioma cell lines (SHG44, A172, U251, and T98G). Lower panels show the same cell-line data re-plotted on an expanded y-axis to visualize differences among glioma cell lines. (B) IRAIN mRNA levels in normal brain tissues versus glioma tissues in the TCGA cohort. (C) Kaplan-Meier curves of overall survival (left) and progression-free survival (right) for patients stratified into high- and low-IRAIN expression groups (log-rank test). (D) Association of IRAIN expression with key clinicopathologic features, including 1p/19q codeletion status (codel vs non-codel), IDH mutation (mutant vs wild-type), MGMT promoter methylation (methylated vs unmethylated), histological subtype (GBM vs LGG), patient age (≤65 vs >65 years), and sex (female vs male). (E) Gene Ontology (left) and KEGG pathway (right) enrichment analyses of IRAIN-correlated genes, highlighting chromatin organization and histone modification-related terms, as well as cancer- and metabolism-related signaling pathways (e.g., FoxO, JAK-STAT, AMPK, mTOR, and cellular senescence pathways). Data are presented as mean ± SD where applicable. *P < 0.05, **P < 0.01, ***P < 0.001. Please click here to view a larger version of this figure.

Figure 3: Identification of immunometabolic-related genes (IMRGs) associated with glioma progression. (A) Volcano plot of differentially expressed genes between high- and low-IRAIN expression groups in the TCGA cohort (|log2FC| > 1.5, FDR < 0.05). Upregulated and downregulated genes are shown in red and blue, respectively; grey dots denote non-significant genes. (B) Determination of the soft-thresholding power for WGCNA. Left: scale-free topology model fit (signed R2) as a function of the soft-threshold power; β = 8 (red point) meets the approximate scale-free criterion (R2 ≈ 0.85). Right: corresponding mean connectivity for each power value. (C) Left: module-trait correlation heatmap showing the relationships between module eigengenes and the IRAIN-associated trait (Control vs Treat/high vs low IRAIN status). The brown module exhibits the strongest positive correlation, whereas the cyan and midnight-blue modules are negatively correlated (numbers in each cell indicate Pearson r and P value). Right: scatter plot of gene significance for the trait versus module membership in the brown module (cor = 0.73, P < 1×10-200). (D) Venn diagram showing the overlap among DEGs, WGCNA brown-module genes, and curated immune/metabolism-related gene sets. The 59 overlapping genes are defined as immunometabolic-related genes (IMRGs). (E) Forest plot of univariate Cox regression analyses for IMRGs in glioma. Squares represent hazard ratios with 95% confidence intervals; genes with P < 0.05 are considered prognostically significant. Please click here to view a larger version of this figure.

Figure 4: Construction and performance of the 17-gene immunometabolic risk signature. (A) Heatmap of concordance indices (C-indices) for 101 prognostic models generated from 10 machine-learning algorithms across the TCGA, CGGA693, CGGA325, and GSE43378 cohorts. The RSF + Elastic Net model with α = 0.3 (top row) shows the highest and most consistent C-indices and was selected as the final model.(B) Regression coefficients of the 17 IMRGs in the final RSF + Elastic Net model; yellow bars indicate risk genes (HR > 1) and blue bars indicate protective genes (HR < 1).(C-F) Kaplan-Meier survival curves (left) and time-dependent ROC curves (right) for the IMRG-based risk score in the TCGA (C), CGGA693 (D), CGGA325 (E), and GSE43378 (F) cohorts. Patients are divided into high- and low-risk groups according to the median risk score in each dataset. Log-rank P < 0.001 for all cohorts. Please click here to view a larger version of this figure.

Figure 5: Independent prognostic value of the IMRG risk score and nomogram construction. (A) Univariate Cox regression analysis of overall survival for clinical variables (grade, sex, age, IDH status, 1p/19q status, MGMT methylation) and the IMRG risk score in the TCGA cohort. (B) Multivariate Cox regression showing that the IMRG risk score remains an independent predictor of overall survival after adjustment for clinicopathologic covariates.(C) Association between the IMRG risk score and clinical/molecular features, including MGMT promoter methylation, age, sex, WHO grade, IDH mutation status, and 1p/19q codeletion status.(D) Nomogram integrating the IMRG risk score with clinical variables to predict 1-, 3- and 5-year overall survival (upper left); ROC curves comparing the discriminative ability of the nomogram, risk score and individual clinical factors (upper right); calibration plots for 1-, 3- and 5-year survival (lower left); and time-dependent C-index curves for the risk score and clinical variables (lower right). The nomogram shows high predictive accuracy (C-index = 0.879, 95% CI 0.855-0.903). Please click here to view a larger version of this figure.

Figure 6: Distribution of the IMRG risk score, survival status, and gene expression in four glioma cohorts. (A-D) For the TCGA (A), CGGA693 (B), CGGA325 (C), and GSE43378 (D) cohorts, the upper heatmaps depict expression patterns of the 17 IMRGs ordered by increasing risk score, the middle panels show the distribution of the calculated risk score and the cutoff separating low- and high-risk groups, and the bottom scatter plots indicate survival status (alive vs dead) of each patient. In all cohorts, deaths are predominantly clustered in the high-risk group, supporting the robustness of the IMRG-based risk model. Please click here to view a larger version of this figure.

Figure 7: Immune microenvironment features associated with the IMRG risk score and IRAIN expression. (A) Boxplots comparing ssGSEA-derived immune cell infiltration and immune-related functional scores between low- and high-risk groups.(B) Boxplots showing differences in the fractions of tumor-infiltrating immune cells between IRAIN-low and IRAIN-high groups.(C) Violin plots of stromal, immune, and ESTIMATE scores in low- versus high-risk groups (left) and IRAIN-low versus IRAIN-high groups (right).(D) LinkET correlation network depicting associations between IRAIN, BIRC5, and multiple immune cell subsets; the color scale indicates correlation coefficients and significance.(E) Expression of immune checkpoint genes (e.g., BTLA, LAG3, HAVCR2, CTLA4, PDCD1, PDCD1LG2, SIGLEC15) in low- and high-risk groups.(F) Violin plot of TIDE scores in low- versus high-risk patients, indicating stronger predicted immune evasion in the high-risk group. P values are shown in the plots; *P < 0.05, **P < 0.01, ***P < 0.001. Please click here to view a larger version of this figure.

Figure 8: Somatic mutation landscape and pathway enrichment analyses of the IMRG risk signature. (A) Oncoplots showing the somatic mutation profiles of key genes in the low-risk (left) and high-risk (right) groups in the TCGA cohort. Stacked bars on the right indicate mutation frequencies for each gene. Lower panels: Kaplan-Meier curves comparing overall survival between TMB-high and TMB-low patients (left) and among four combined categories (high-TMB/high-risk, high-TMB/low-risk, low-TMB/high-risk, low-TMB/low-risk; right). (B) GeneMANIA interaction network of the 17 IMRGs and related genes. Node size reflects connection degree; edges are colored according to interaction type (physical interaction, co-expression, co-localization, shared pathways, or predicted interactions). (C) Gene Ontology (left) and KEGG pathway (right) enrichment analyses of genes within the IMRG-related network, highlighting terms associated with antigen processing and presentation, regulation of NF-κB signaling, glutathione metabolism, viral infection pathways, and JAK-STAT signaling. (D) STRING protein-protein interaction (PPI) network of IMRGs (left) and a simplified subnetwork emphasizing hub genes (right), including EGFR, BIRC5, RRM1, IDH1, BMPR1A, CDKN2B, RHOA, PDIA3, CALR, and PRIM1. Please click here to view a larger version of this figure.

Figure 9: IRAIN suppresses glioma proliferation and angiogenesis via the IGF1R-JAK2-STAT3-BIRC5 axis. (A) Cell-viability curves (CCK-8/MTT) for control, IRAIN, IRAIN + si-IGF1, and IRAIN + SURVIVIN groups measured at 0, 24, 48, and 72 h. (B) Representative Annexin V-FITC/PI flow-cytometry plots showing apoptosis in the indicated groups. (C) Immunocytochemical staining for CD31 and quantification of the immunoreactive score, illustrating reduced endothelial marker expression in IRAIN-overexpressing and IRAIN + si-IGF1 cells, with recovery after SURVIVIN overexpression. (D) Western blot analysis and densitometric quantification of IGF1, IGF1R, JAK2, STAT3, and SURVIVIN/BIRC5 protein levels in control, IRAIN, IRAIN + si-IGF1, and IRAIN + STAT3 groups. Data are presented as mean ± SD; *P < 0.05, **P < 0.01, ***P < 0.001. Please click here to view a larger version of this figure.

Table 1: Primer sequences used for qRT-PCR. Melt-curve analyses confirmed single specific products under the qRT-PCR conditions described in the protocol.

Supplementary File 1: Raw data.zip. Please click here to download this File.

Glioma remains one of the most lethal malignancies of the central nervous system, characterized by profound intratumoral heterogeneity and resistance to conventional therapy. Despite surgical resection combined with chemoradiotherapy, recurrence and mortality remain high, and the median survival of patients, particularly those with glioblastoma, has shown limited improvement over recent decades^4,24. Increasing evidence indicates that metabolic reprogramming and immune suppression within the tumor microenvironment (TME) are key determinants of glioma aggressiveness^15,25. Yet, patients with similar clinical characteristics often display substantial molecular and metabolic heterogeneity, underscoring the need for integrated biomarkers that capture both immune and metabolic dimensions of glioma biology.

In this study, we identified and functionally characterized immunometabolic-related genes (IMRGs) associated with IRAIN dysregulation. Using differential expression, WGCNA, and Cox regression analyses, we established a 17-gene IMRG prognostic signature through an integrated leave-one-out cross-validation framework combining ten machine-learning algorithms and 101 model combinations. The optimal Random Survival Forest-Elastic Net (α = 0.3) model exhibited strong predictive accuracy across TCGA, CGGA, and GEO cohorts, outperforming established clinicopathologic indicators such as WHO grade, IDH mutation, MGMT methylation, and 1p/19q codeletion. Further GO/KEGG enrichment analysis revealed that the screened immune and metabolic pathways co-expressed with IRAIN in the JAK-STAT pathway. IRAIN downregulation inhibited IGF1R activity in glioma. Then, IGF1R can engage JAK kinases to activate STAT3, a transcription factor that coordinates tumor survival, immunosuppression, and immunometabolic programs in glioma^20,21. Consistent with this biology, several canonical STAT3 effectors represent plausible points of convergence between IRAIN loss and components of the IMRGs landscape. Furthermore, experimental results proved that IRAIN overexpression suppressed the activity of the IGF1R-JAK2-STAT3 pathway to reduce BIRC5 expression, leading to a decrease in tumor vessel formation, apoptosis, and angiogenesis in glioma stem cells. These findings establish a mechanistic bridge between transcriptomic reprogramming and phenotypic tumor suppression, confirming that IRAIN acts as a key upstream regulator within the immunometabolic landscape of glioma. Importantly, this integrated framework links a large-scale bioinformatics signature to experimentally validated molecular mechanisms. Different regions of the TME exhibit distinct metabolic characteristics. Metabolic alterations can induce antitumor immune suppression by modulating various metabolic pathways in immune cells, thereby influencing tumor heterogeneity and survival duration^26,27. Tumor stem cells appear to exploit this altered environment, facilitating their physiological invasion of normal glioma regions²⁸. Metabolic adaptation in glioma is intimately intertwined with immune dysfunction¹⁵. Tumor-derived metabolites such as lactate, arginine, and glutamine reshape immune cell function and cytokine production, promoting tolerance and T-cell exhaustion^29,30. The hypoxic and acidic microenvironment fosters an immunosuppressive niche dominated by Tregs and M2 macrophages^31,32,33, as reflected by our observation that low IRAIN expression correlates with elevated immune and stromal scores, high TIDE values, and increased infiltration of immunosuppressive cell types. This metabolic-immune convergence explains why the IMRG high-risk subgroup exhibits both enhanced metabolic activity and impaired antitumor immunity. The integration of metabolic and immune signaling, particularly through the JAK-STAT-BIRC5 axis, suggests that IRAIN loss may contribute to immune evasion and therapy resistance by reinforcing survival and angiogenic signaling pathways.

The clinical implications of these findings are twofold. First, the 17-gene IMRG signature provides a robust framework for prognostic stratification and therapeutic guidance in glioma. Its superior predictive power relative to traditional histopathologic features underscores the potential of incorporating immunometabolic indices into future clinical decision algorithms. Second, IRAIN emerges as a mechanistic regulator linking transcriptional, metabolic, and immune pathways. By repressing IGF1R transcription and downstream STAT3 activation, IRAIN inhibits the expression of BIRC5, a potent anti-apoptotic protein that mediates resistance to apoptosis, angiogenesis, and chemotherapeutic stress^34,35,36. These data align with prior studies in breast, laryngeal, and lung cancers showing that antisense lncRNA IRAIN suppresses IGF1R expression through promoter-enhancer interactions, thereby attenuating oncogenic signaling^17,18,19,37. Furthermore, our work extends this mechanism to glioma, demonstrating that IRAIN deficiency promotes survival signaling via JAK-STAT-BIRC5 activation. Western blotting and functional assays confirmed that IRAIN overexpression reduces p-JAK2 and p-STAT3 levels, downregulates BIRC5, and enhances apoptosis, while BIRC5 overexpression reverses these effects. Given that IGF1R inhibitors and JAK/STAT blockers (e.g., AG490) are already under clinical evaluation³⁸, the IRAIN-IGF1R-JAK-STAT-BIRC5 axis represents a tractable therapeutic target for future combinatorial strategies. Collectively, these results highlight the convergence of immune and metabolic dysregulation as a hallmark of glioma progression and identify IRAIN as a central node in this network. The integrated IMRG model, validated across multiple datasets, not only provides a reliable prognostic biomarker but also offers mechanistic insights into how metabolic reprogramming and immune evasion are co-regulated. Targeting this regulatory axis may enhance immunotherapy responsiveness and improve clinical outcomes in glioma patients.

In conclusion, our study defines IRAIN as a master regulator of glioma immunometabolism that suppresses tumor growth by modulating the IGF1R-JAK2-STAT3-BIRC5 pathway. The 17-gene IMRG signature serves as a clinically applicable tool for risk stratification, while mechanistic validation confirms the biological relevance of this network. Together, these findings bridge computational discovery with functional biology and open avenues for developing metabolism- and immunity-oriented therapeutic interventions in glioma.