Clinical information and transcriptomic data of LUAD patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), and data from the GES31210 dataset, including transcriptomic, phenotypic, and clinical details, were accessed from the Gene Expression Omnibus (GEO) (NCBI, http://www.ncbi.nlm.nih.gov/geo/). Data were processed in R software (version 4.1.3) using the “limma” package for normalization and log transformation. GREM1 expression levels in normal and tumor tissues were compared using the Wilcoxon test, with statistical significance set at two-tailed P < 0.050. To control for potential false positive results arising from testing a large number of genes simultaneously, we applied the Benjamini-Hochberg False Discovery Rate (FDR) correction. This approach ensured the statistical reliability of the identified differentially expressed genes while minimizing the false positive rate to the greatest extent possible.

We conducted a WGCNA on preprocessed gene expression data using the “WGCNA” package in R software (version 4.1.3). Outlier samples were identified and removed as an initial step. Subsequently, we employed the pickSoftThreshold function to test and select an appropriate soft threshold. Gene co-expression modules were identified using the dynamic tree cut method, and the module eigengenes (MEs) for each module were calculated. Finally, Pearson correlation analysis and Student’s t-tests were conducted to evaluate relationships between module eigengenes and patient clinical characteristics.

Gene annotation was performed using the “org.Hs.eg.db” package in R, with the additional installation of “ggplot2”, “colorspace”, “stringi”, “clusterProfiler”, and “enrichplot” packages. The “enrichGO” function from the “clusterProfiler” package was used for GO enrichment analysis, while the “enrichKEGG” function was used for KEGG enrichment analysis. In the GO enrichment analysis and KEGG pathway functional enrichment analysis, we tested the significance of multiple biological pathways and functions. To address statistical errors caused by multiple hypothesis testing, we also applied FDR correction to adjust the p-values. This method ensured the statistical robustness of the enrichment analysis results. The results of the enrichment analyses were visualized using circular plots.

Patients were classified into high- and low-risk groups according to the median expression level of the GREM1 gene. The association between GREM1 expression and OS was illustrated with Kaplan-Meier survival curves, and the log-rank test was applied to calculate survival differences between the two groups. Univariate and multivariate Cox regression analyses were performed using the “survival” and “forestplot” R packages to calculate hazard ratios (HR) and 95% confidence interval (CI), enabling assessment of GREM1’s prognostic value. The Cox regression model included GREM1 expression levels and other clinical covariates potentially influencing survival (e.g., age, gender, and tumor stage) to assess its independent prognostic value through multivariate analysis. Since survival analysis primarily focused on single hypothesis testing, multiple comparison corrections were not applied.

Based on the expression levels of GREM1 and clinical stage, a nomogram was constructed using the “rms” R package to predict the OS probabilities for patients with LUAD at 1, 3, and 5 years. This method provides an effective and convenient approach for estimating individual survival rates in patients. Calibration curves were also plotted to assess the nomogram’s predictive accuracy.

The “estimate” R package and the ESTIMATE algorithm were used to analyze differences in stromal and immune scores across samples. Bioinformatics tools, including TIMER 2.0 and CIBERSORT, assessed the levels of 22 infiltrating immune cell subtypes. The “limma” and “CIBERSORT” R packages were used to quantify specific cell types in the samples based on overall expression data. In exploring the relationship between GREM1 gene expression and immune cell infiltration (such as T cells and macrophages) as well as immune scores (e.g., ESTIMATE scores), Spearman correlation coefficients were calculated. To minimize the impact of multiple comparisons on the results, P-values from the correlation analysis were adjusted using FDR correction.

We used the “maftools” package in R to process single-nucleotide variant (SNV) data from the TCGA database and calculated each sample’s TMB and MSI. Spearman correlation analysis was conducted to determine the association between GREM1, TMB, and MSI. Additionally, the “fmsb” package was used to create radar charts to visualize the correlation analysis results for TMB and MSI.

Fresh tissue samples were collected from five LUAD patients who underwent surgical treatment at Harbin Medical University Cancer Hospital. None of the patients received anticancer therapy prior to surgery. This study was approved by the Ethics Committee Review Board of the Harbin Medical University Cancer Hospital and in conformity to the Declaration of Helsinki. Informed written consent was obtained from all patients before participation in the study.

Human LUAD cell lines (PC-9, A549, H1299, H1650, HCC827) and HBE cells were obtained from the Laboratory of Medical Genetics (Department of Biology, Harbin Medical University, Harbin, China). The HBE, A549, H1299, H1650, and HCC827 cell lines were cultured in RPMI 1640 (Gibco, Life Technologies, California, USA), while the PC-9 cell line was cultured in DMEM (Gibco, Life Technologies, California, USA). Both culture media were supplemented with 10% fetal bovine serum (ScienCell, California, USA) and 1% penicillin-streptomycin (PYG0016, Boster). Cell culture was performed in a humidified incubator at 37°C with 5% CO2.

Total protein was extracted from LUAD specimens and cells using RIPA buffer (Beyotime, China). Equal amounts of protein were subjected to electrophoresis on 10% SDS-PAGE and transferred to a PVDF membrane. The membrane was blocked with 5% non-fat dry milk (Beyotime, China) in TBST for 1 hour at room temperature, followed by overnight incubation at 4°C with primary antibodies against GREM1 (1:1000, Affinity), E-cadherin (1:20000, Proteintech), N-cadherin (1:3000, Proteintech), Vimentin (1:50000, Proteintech), and GAPDH (1:1000, Absin). After washing with TBST, the membrane was incubated with secondary antibodies (1:10000, Beijing Zhongshan Golden Bridge Biotechnology Co., Ltd.) for 1 hour at room temperature. Finally, target protein bands were detected using the ECL detection system (Beyotime, China), with GAPDH serving as the internal controls.

A549 cells were transfected with GREM1-OE (RiboBio, China) according to the manufacturer’s instructions to upregulate GREM1 expression. H1650 cells were transfected with GREM1 siRNA (RiboBio, China) as per the manufacturer’s guidelines to downregulate GREM1 expression.

The proliferation of A549 and H1650 cells was assessed using the Cell Counting Kit-8 (CCK-8) (APExBIO, USA). Approximately 4 × 10³ cells were seeded into 96-well plates and cultured for 24, 48, 72, and 96 hours. Optical density (OD) at 450 nm was measured with a microplate reader.

In the wound healing assay, when the cells reached confluence, we created wounds in the cell layer using a 1000 µL pipette tip. Serum-free medium was then applied to maintain cell viability. Wound areas were photographed immediately after wounding (0 hours) and 24 hours later using an optical microscope (Nikon, Japan).

In the Transwell assay, 2 × 10⁴ cells were seeded into each upper chamber of 24-well plates, with or without Matrigel (Becton Dickinson). After 24 or 48 hours of incubation (with or without Matrigel), cells on the bottom membrane were fixed with paraformaldehyde, stained with crystal violet, and analyzed microscopically.

Statistical analyses were conducted using R version 4.1.3. Survival curves were generated with the Kaplan-Meier method and analyzed statistically. The Wilcoxon test was used to compare differences between groups, while Pearson and Spearman correlation coefficients were calculated for association analyses. To control for multiple comparisons, the Benjamini-Hochberg False Discovery Rate (FDR) correction was applied when analyzing differentially expressed genes and pathway enrichment analyses. Statistical significance was defined as a two-tailed P value of < 0.050.

First, we summarized the research design workflow, including the screening process of GREM1 and the exploration of its biological function ( Figure 1 ).

Differential gene expression analysis was performed between LUAD tissues and normal tissues using the TCGA and GSE31210 datasets. In the TCGA dataset, 4,859 differentially expressed genes were identified (|log2(FC)| > 1, adj.P.Val < 0.05), of which 2,347 were upregulated and 2,512 were downregulated. Similarly, in the GSE31210 dataset, 3,315 differentially expressed genes were detected (|log2(FC)| > 1, adj.P.Val < 0.05), with 1,785 upregulated and 1,530 downregulated. To identify key shared differentially expressed genes across both datasets, we conducted an intersection analysis, revealing 1,927 common genes. We then integrated clinical data, including age, sex, and disease stage, and performed weighted gene co-expression network analysis (WGCNA) on 1,927 differentially expressed genes. The scale-free fit index and average connectivity were calculated to determine a soft threshold of β=4 ( Figures 2A, B ). Using the dynamic tree cut method, we identified seven gene co-expression modules, as shown in the clustering tree in Figure 2C . The module-trait correlation heatmap ( Figure 2D ) revealed a significant correlation between the green module and the clinical stage of LUAD (cor=0.33, P=5e-07), suggesting that the 72 genes in this module may be closely associated with LUAD invasion and metastasis. Consequently, the green module was identified as the key module for further investigation.

To further explore the biological functions of genes in the green module, we conducted Gene Ontology (GO) analysis ( Figure 2E ). The results of the biological process (BP) analysis revealed that these genes are strongly associated with the organization of the extracellular matrix and extracellular structures. In terms of cellular components (CC), these genes participate in the formation of collagen-containing extracellular matrices and the lumen of the endoplasmic reticulum. Regarding molecular functions (MF), these genes contribute to activities such as the structural composition of extracellular matrix components. Through KEGG enrichment analysis ( Figure 2F ), we observed that genes within the green module are significantly enriched in essential biological pathways, such as protein digestion and absorption, ECM-receptor interaction, and the PI3K-Akt signaling pathway. These findings indicate that genes in the green module may impact tumor cell invasion and metastasis by modulating the formation and function of the TME, thus playing a crucial role in the progression and pathology of LUAD.

In this study, we selected the top 25 genes with the highest connectivity from the 72 genes in the green module and performed an intersection analysis with 2,660 genes from the immune database. This analysis identified two candidate genes: GREM1 and PLAU. To evaluate the prognostic impact of these two genes on LUAD patients, we conducted univariate Cox regression analysis ( Figure 2G ). The results showed that, of the two candidate genes, only the expression level of GREM1 was significantly associated with LUAD patient prognosis. (P = 0.020, HR = 1.126, 95% CI: 1.019–1.245).

We assessed the prognostic impact of GREM1 in LUAD patients across multiple datasets. In the TCGA-LUAD dataset, which includes 594 patient samples (535 tumor samples and 59 normal samples), GREM1 expression levels were significantly higher in tumor tissues compared to normal tissues (P<0.001). Similarly, in the LUAD dataset GSE31210, which comprised 246 patient samples (226 tumor samples and 20 normal samples), GREM1 expression was significantly higher in tumor tissues (P<0.001) ( Figure 3A ). Additionally, we conducted Kaplan-Meier survival analysis to compare overall survival (OS) between patients with high and low GREM1 expression in the TCGA and GSE31210 datasets. The analysis revealed that patients with high GREM1 expression had significantly shorter OS (P<0.05) ( Figure 3B ).

To further assess the prognostic and diagnostic significance of GREM1 in LUAD, we conducted univariate and multivariate Cox regression analyses. Univariate Cox analysis revealed a significant association between high GREM1 expression and prognosis. In the multivariate Cox proportional hazards regression model for the TCGA and GSE31210 datasets, high GREM1 expression was confirmed as an independent prognostic biomarker (for OS, HR=1.803 and 2.286, 95% CI=1.185-2.741 and 1.569-3.330, P=0.006 and <0.001, respectively; for disease-specific survival (DSS), HR=1.988, 95% CI=1.142-3.461, P=0.015) ( Figure 3C ).

To improve prediction of the clinical prognostic value of GREM1, we incorporated GREM1 expression and tumor staging-both independent mortality predictors-into a multivariable Cox proportional hazards model. Using this model, we developed a prognostic nomogram to predict the OS of LUAD patients at 1, 3, and 5 years ( Figure 3D ), achieving predictive accuracies of 0.897, 0.618, and 0.364, respectively. This model visually represents the combined impact of GREM1 expression and tumor staging on LUAD prognosis. We evaluated and validated this nomogram using TCGA data. Calibration curve analysis confirmed that the model provides reliable accuracy in predicting 1-, 3-, and 5-year OS, aligning with the standards of an ideal predictive model ( Figure 3E ).

Initially, we evaluated the role of GREM1 in the invasion and metastasis of LUAD. We categorized the samples using the MEXPRESS tool based on several clinical factors, including the occurrence of new tumors post-treatment, tumor stage, smoking history, and sample type. The results revealed significant differential expression of GREM1 related to tumor status (P=0.03) and staging (P=0.03), with significantly higher expression levels in stages T2 and T3 ( Figure 4A ). Furthermore, GREM1 expression was significantly higher in stage II and III patients than in stage I patients ( Figure 4B ). This finding is consistent with our previous studies, which suggest that GREM1 expression is closely linked to LUAD tumor progression. Although the GREM1 expression differences in T4 and stage IV did not reach statistical significance, this may be attributed to an insufficient sample size.

We next analyzed GREM1 expression across six immune subtypes to investigate its relationship with the intratumoral immune state ( Figure 4C ). The results revealed a significant upregulation of GREM1 expression in the wound healing, IFN-gamma dominant, and immune-related inflammation subtypes. These findings suggest that GREM1 is not only involved in the invasion and metastasis of LUAD but also closely associated with the TME and immune response.

To further explore the impact of GREM1 expression on the TME, we compared its correlation with the ESTIMATE score, Immune score, and Stromal score. GREM1 expression in the TCGA and GSE31210 datasets was significantly positively correlated with the immune score (P=3.6e-07, R=0.24; P=0.00037, R=0.23), stromal score (P<2.2e-16, R=0.6; P=1.2e-12, R=0.45), and estimate score (P<2.2e-16, R=0.44; P=3.4e-09, R=0.38) ( Figures 4D, E ). Additionally, significant differences in these scores were observed between the high and low GREM1 expression groups ( Figures 4F, G ). These results suggest that GREM1 overexpression may promote tumor growth and metastasis by modifying the physical or chemical properties of the TME.

A significant positive correlation was observed between GREM1 mRNA expression levels and Log2 copy number variation, suggesting that increased copy numbers are associated with elevated GREM1 expression (P=0.0235, R=0.10) ( Figure 4H ). Further analysis of GREM1 gene alterations in TCGA-LUAD tissues, including deep and shallow deletions, diploidy, and copy number gains, indicated a marked increase in GREM1 copy number ( Figure 4I ).

Based on the above results, we applied the MCP counter algorithm to examine the association between GREM1 and immune cells. Comparative analysis of the TCGA and GSE31210 datasets demonstrated a significant association between GREM1 expression and the infiltration of multiple immune cells ( Figure 5A ). Following this, the relationship between GREM1 expression and 22 immune cell types was analyzed using the CIBERSORT algorithm ( Figure 5B ). Notably, significant differences in the infiltration levels of specific immune cells, including CD4 memory resting T cells, CD4 memory activated T cells, activated NK cells, monocytes, macrophages, resting mast cells, and neutrophils, were observed between high and low GREM1 expression groups. These findings indicate that GREM1 may affect immune cell composition within the tumor immune microenvironment (TIME).

To explore the potential gene networks or pathways associated with GREM1, we examined its correlations with other genes. As shown in Figure 5C , in the TCGA dataset, GREM1 exhibited a significant positive correlation with twist family BHLH transcription factor 1 (TWIST1), actin alpha cardiac muscle 2 (ACTA2), and platelet-derived growth factor receptor alpha (PDGFRA) (R>0.6, P<0.001), and a significant negative correlation with SMAD family member 9 (SMAD9) (R<-0.2, P<0.001). In the GSE31210 dataset, GREM1 showed a significant positive correlation with TWIST1 (R>0.6, P<0.001), and a significant negative correlation with transforming growth factor beta receptor 2 (TGFBR2), T-Box transcription factor 2 (TBX2), and SMAD9 (R<-0.5, P<0.001). The chemokine correlation analysis ( Figure 5D ) showed that GREM1 had a strong positive correlation with CXCL12, CCL8, and CCL21 in the TCGA dataset (R>0.4, P<0.001). Similarly, in the GSE31210 dataset, GREM1 was positively correlated with CXCL10 and CXCL11 (R>0.4, P<0.001).

Furthermore, we examined the effect of GREM1 on clinical responses to ICIs in LUAD. Correlation analysis revealed that higher GREM1 expression was strongly associated with elevated immune function scores, such as C-C chemokine receptor (CCR), checkpoint markers, inflammation-promoting factors, macrophages, T-cell co-inhibition, tumor-infiltrating lymphocytes (TILs), and regulatory T cells (Tregs) ( Figure 6A ). Furthermore, patients with high GREM1 expression exhibited a significantly higher tumor mutation burden (TMB) than those with low expression (P=1.2e-06) ( Figure 6B ). Additionally, as GREM1 expression increased, the RNA stemness score (RNAss) gradually decreased (P=0.011, R=-0.12) ( Figure 6C ). In conclusion, GREM1 may play a critical role in immune cell recruitment and the response to immunotherapy.

Considering the current scarcity of research on GREM1 in cancer studies, particularly regarding solid tumors, we conducted a systematic evaluation of GREM1 expression across various solid tumor types. Results indicate that GREM1 is highly expressed in LUAD and is overexpressed in other cancers such as BRCA, HNSC, KIRC, KIRP, THCA, and UCEC ( Figure 6D ). Univariate Cox analysis further confirmed the prognostic significance of GREM1 in multiple solid tumors beyond LUAD, such as KIRC, KIRP, MESO, PAAD, and STAD ( Figure 6E ). Analysis of GREM1’s TMB and microsatellite instability (MSI) revealed a strong association between its expression and both the proliferative and metastatic potential of LUAD, along with its response to immunotherapy ( Figures 6F, G ).

We conducted in vitro experiments to investigate the pathogenic mechanism of GREM1 in LUAD. First, we assessed GREM1 expression in tissue samples using Western blot analysis. The results showed significantly higher GREM1 expression in cancerous tissues compared to adjacent non-cancerous tissues in five paired samples ( Figure 7A ). We then collected various cell lines, including LUAD cell lines (PC-9, A549, H1299, H1650, HCC827) and a human bronchial epithelial cell line (HBE), to analyze GREM1 expression levels using Western blot. The results demonstrated generally higher GREM1 expression in LUAD cell lines compared to normal cells ( Figure 7B ). Based on these findings, we selected the H1650 cell line, which showed the highest GREM1 expression, and the A549 cell line, which showed the lowest expression, for further experiments.

A549 cells were transfected with the GREM1 overexpression plasmid. Western blot analysis confirmed that GREM1 levels were significantly higher in the GREM1 overexpression group compared to the control group. A GREM1-targeting siRNA was designed and transfected into H1650 cells, while non-targeting siRNA was used as a control. Western blot results showed that GREM1 expression in the knockdown group was significantly lower than in the control group ( Figure 7C ).

Figure 7D shows the effect of GREM1 on EMT-related proteins in A549 and H1650 cells as observed through Western blot analysis. The effects of GREM1 overexpression or silencing on the expression of the epithelial marker E-cadherin and mesenchymal markers N-cadherin and Vimentin were examined. In A549 cells, overexpression of GREM1 significantly downregulated E-cadherin while upregulating N-cadherin and Vimentin, suggesting that GREM1 overexpression promotes the EMT process. Conversely, in H1650 cells, GREM1 knockdown significantly upregulated E-cadherin and downregulated N-cadherin and Vimentin, indicating that GREM1 inhibition may hinder EMT. These findings highlight GREM1’s regulatory role in EMT: its overexpression promotes EMT, while its knockdown suppresses the process. This suggests that GREM1 may influence tumor cell migration and invasion through EMT modulation.

CCK-8 assays were conducted on A549 and H1650 cell lines to evaluate the effect of GREM1 on cell proliferation. The results showed that GREM1 overexpression significantly promoted cell proliferation in A549 cells compared to the control group. Conversely, GREM1 knockdown in H1650 cells significantly decreased proliferation (P<0.01) ( Figure 7E ), suggesting that GREM1 is crucial for tumor cell proliferation.

Wound healing assay results indicate that GREM1 overexpression significantly promotes cell migration in A549 cells. At 24 hours, the wound closure rate in the GREM1 overexpression group was significantly higher than in the control group, suggesting that GREM1 enhances A549 cell migration. Similarly, in H1650 cells, GREM1 knockdown significantly reduced cell migration, with the wound closure rate in the siGREM1 group lower than in the control group at 24 hours, indicating that GREM1 downregulation suppresses H1650 cell migration ( Figure 7F ).

Transwell assay results demonstrated that GREM1 expression in A549 and H1650 cells is closely associated with invasion and migration. In A549 cells, GREM1 overexpression significantly enhanced invasive and migratory abilities, as evidenced by more cells passing through the chamber in the overexpression group than in the control group. Conversely, in H1650 cells, GREM1 knockdown significantly inhibited invasion and migration, with fewer cells passing through the chamber in the siGREM1 group than in the control group ( Figure 7G ). These findings suggest that GREM1 positively regulates cell invasion and migration.