Samples of lung adenocarcinoma tissue (n = 75) and normal lung tissue (n = 50) were obtained from patients at the Fujian Medical University Union Hospital between September 2018 and December 2019. All samples were confirmed by pathological diagnosis and none of the donor patients had received any other type of treatment prior to surgery. The study protocol was approved by the Ethics Committee of Fujian Medical University Union Hospital, and all patients provided their written informed consent.

Four human LUAD cell lines (H292, H1975, H1299, and A549) and a normal lung epithelial cell line (BEAS-2B) were acquired from the Cell Bank of the Chinese Academy of Sciences. All the cell lines were cultured in RPMI-1640 medium (Gibco, Waltham, MA, USA) containing 10% fetal bovine serum (FBS, Gibco) and 1% penicillin-streptomycin (HyClone, Logan, UT, USA) in a humidified 37°C incubator with a 5% CO₂ atmosphere.

Small interfering RNA sequences against GCC2-AS1 (si-GCC2-AS1) and a negative control (NC) sequence were obtained from GenePharma (Shanghai, China). For siRNA transfection, cells were cultured and then transfected with si-GCC2-AS1 or the NC using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) as described in the manufacturer’s instructions.

Total RNA was extracted using TRIzol reagent (Invitrogen) and then reverse transcribed to cDNA using a PrimeScript™ RT reagent kit (Takara, Shiga, Japan) according to manufacturer’s instructions. GCC2-AS1 expression was determined by using SYBR Green GoTaq Master Mix (Promega, Madison, WI, USA) on an ABI 7500 thermocycler (Applied Biosystems Foster City, CA, USA). The specific primers used for GCC2-AS1 (Forward: ACCTGGTCTGGATCGGTCAC; Reverse: GGGCATGTTCTTCTTGACTGC) were synthesized by GenePharma (Shanghai, China). An 18S ribosomal RNA was used as internal control, and 2 ^{-ΔΔ Ct} method was utilized to calculate the relative expression of GCC2-AS1 (14).

Digoxin (DIG)-labeled GCC2-AS1, a modified probe for FISH, was synthesized by Servicebio (Wuhan, China). The sequence of the probe was 5’-DIG-CCGCCATCCTTGTTGTAGGGAACTCG-DIG-3’. A549 and H1299 cells were seeded into a confocal dish, and the subsequent procedures were performed as previously described (15). Images were captured with a confocal microscope (Nikon, Japan).

Cells were seeded into the wells of 96-well plates (3,000 cells/well) containing 200 μl of RPMI-1640 medium. The medium was refreshed every 24 h. Next, CCK8 solution (10 μl) was added to each well, and the cells were incubated for 2 h. Finally, the absorbance of each well at a wavelength of 450 nm was measured with a microplate reader.

Cell invasion assays were performed using a 24-well Transwell chamber (Corning, NY, USA) coated with Matrigel (Corning). Cells (1.5 × 10⁵) suspended in 1% FBS medium were added to the upper chamber and 600 μl of 10% FBS medium was added to the lower chamber. After 24 h of incubation, cells that had invaded the lower chamber were fixed with 4% paraformaldehyde (PFA) and stained with 1% crystal violet.

RNA-seq data and corresponding clinical information obtained from various LUAD projects (513 patients, Workflow Type: HTSeq-FPKM) were downloaded from the TCGA database and subjected to further analysis. Later, level 3 HTSeq-FPKM data was transformed into TPM (transcripts per million reads) data for subsequent analysis. Some cases with unavailable or unknown clinical characteristics were considered as missing values. Data containing a differential expression analysis of GTEx (Genotype-Tissue Expression) or a pan-cancer analysis were retrieved from University of California Santa Cruz (UCSC) XENA (https://xenabrowser.net/datapages/). RNAseq data for TPM obtained from the TCGA and GTEx analysis were continuously processed using Toil software (16). Based on GCC2-AS1 expression, the median expression value was used as a cut-off value, and tumor samples were then categorized into low and high groups, respectively. This study was conducted in accordance with publication guidelines provided by TCGA (https://cancergenome.nih,gov/publications/publicationguidelines).

The DESeq2 R package was used to obtain expression profiles (HTSeq-Counts) that were subsequently used to make comparisons between the low and high GCC2-AS1 expression groups. A log (fold change) ||log (FC)|>1.5 and an adjusted P-value<0.05 were considered to constitute a threshold value for DEGs. Volcano plots and a heat map of DEGs were created using the gglot2 package (version 3.1.0, https://gihub.com/tidyverse/ggplot2) and pheatmap package (version 1.0.10, https://CRAN.R-project.org/package=pheatmap), respectively.

The clusterProfiler package (17) was used to perform a gene ontology (GO) enrichment analysis of the input gene list, which included Biological Process (BP), Cellular Components (CC), and Molecular Function (MF). The same package was used to perform a KEGG pathway enrichment analysis. Furthermore, a gene set enrichment analysis (GSEA) was also performed using the clusterProfiler R package to determine the critical functions and pathways that involved the previously defined gene sets in the low and high GCC2-AS1 expression groups. Gene set permutations were carried out 1,000 times for each analysis. Based on |Normalized Enrichment Score||NES|>1, P-value <0.05, and false discovery rate (FDR) q value <0.25, a notably enriched pathway was identified for each phenotype.

All experimental data were analyzed using IBM SPSS Statistics for Windows, Version 22.0 software (IBM Corporation, Armonk, NY, USA). Statistical analyses of bioinformatics data were performed using R (v.3.5.1). Wilcoxon rank sum tests were used to compare the levels of GCC2-AS1 expression between lung adenocarcinoma and normal groups. Pearson’s chi-square test was used to investigate the relationship between GCC2-AS1 expression and clinicopathologic parameters. Additionally, a univariate Cox regression analysis and the Kaplan-Meier method were used to analyze the association between clinicopathologic parameters and OS time. We also conducted multivariate Cox regression analyses to identify independent prognostic factors. The individual hazard ratio (HR) and 95% confidence interval (95% CI) were determined to assess the hazard risk of each parameter. Subsequently, valuable clinicopathological parameters identified from multivariate Cox regression analyses were used to construct a nomogram for predicting the prognosis of LUAD patients. The nomogram was constructed using the rms package in R (Version: 5.1-3; https://cran.r-project.org/web/packages/rms/index.html). Additionally, calibration curves were created to check the nomogram’s predictive efficiency. A ROC (receiver operating characteristic) analysis was carried out using the pROC package (18) in R to estimate the model’s ability to discriminate between lung adenocarcinoma and normal groups. All statistical tests were 2-sided and a P-value <0.05 was considered to be statistically significant.

We initially downloaded RNA-seq data from the UCSC XENA databases and then used Toil workflow (16) software to transform it into TPM data for the purpose of detecting the patterns of GCC2-AS1 expression in 33 pan-cancer and normal tissue sample by performing Wilcoxon rank sum tests. As shown in Figures 1A and 1D , levels of GCC2-AS1 expression were compared with normal tissues and adjacent tissues in online database, respectively. Surprisingly, significantly higher levels of GCC2-AS1 expression were found in the LUAD tissues when compared with the normal tissues ( Figure 1B ). An ROC analysis showed that GCC2-AS1 had a high ability to discriminate between LUAD patients and healthy individuals, with an AUC of 0.785 ( Figure 1C ), suggesting that GCC2-AS1 could be used to help identify people with LUAD. Furthermore, we demonstrated that GCC2-AS1 expression was upregulated in the LUAD samples (n = 513) when compared with its expression in samples of adjacent tissue (n = 59, Figure 1E ). We next used the Wilcoxon signed rank test to compare the levels of GCC2-AS1 expression in 57 samples of LUAD tissue with those in paired adjacent tumors. The results showed that GCC2-AS1 was overexpressed in the LUAD samples ( Figure 1F ). These data suggested that GCC2-AS1 might function as an oncogene that facilitates LUAD carcinogenesis.

To determine whether GCC2-AS1 was differentially expressed in LUAD tissues, we collected 75 samples of LUAD tissue and 50 samples of normal lung tissue and analyzed them by RT-qPCR. We found that the levels of GCC2 expression in the LUAD tissues were extremely upregulated when compared those in normal lung tissues ( Figure 2A ) which is consistent with online database we discovered. We also found that GCC2-AS1 was more highly expressed in the LUAD tissues obtained from metastatic patients (n = 23) than from non-metastatic patients (n = 52, Figure 2B ). Moreover, a similar result was found in a comparison of the LUAD cell lines with the normal lung epithelial cell line ( Figure 2C ). Finally, we found that GCC2-AS1 was mainly distributed in the cytoplasm of both A549 and H292 cells ( Figure 2D ).

The clinical characteristics of 513 LUAD patients (276 females and 237 males, mean age = 63 years) and their levels of GCC2-AS1 expression were obtained from the TCGA cohort ( Table 1 ) We then divided those patients into two groups (a high and low expression group GCC2-AS1, respectively) based on the median value for GCC2-AS1 expression. An analysis of the total cohort showed that the parameters of age (P = 0.151), T stage (P = 0.111), N stage (P = 0.577), M stage (P = 0.095), pathologic stage (P = 0.084), primary therapy outcome (P = 0.962), race (P = 0.267), residual tumor (P = 0.721), anatomic neoplasm subdivision (Left/Right, P = 0.272; central/peripheral, P = 0.755), and tumor status (P = 0.311) were not associated with GCC2-AS1 expression. In contrast, a correlation analysis conducted using the chi square test showed that a high level of GCC2-AS1 expression was significantly correlated with gender (P = 0.007), smoker (P = 0.012), Tumor Protein P53 (TP53) mutation status (P<0.001), and Kirsten Rat Sarcoma Viral Oncogene (KRAS) mutation status (P = 0.011). Another analysis conducted using a similar method revealed that a high level of GCC2-AS1 expression was significantly associated with N stage (P = 0.031) and TNM stage (P = 0.005) ( Table 2 ).

To explore the biological effects of GCC2-AS1 on LUAD cells, we decreased the GCC2-AS1 levels in LUAD cells by transient transfection with siRNA. Transfection efficiency was determined by RT-qPCR after 48 h ( Figure 3A ). CCK-8 assays were conducted to detect the influence of GCC2-AS1 on LUAD cell proliferation. The results suggested that decreased GCC2-AS1 levels significantly suppressed the proliferation of LUAD cells when compared with cells in the NC groups ( Figure 3B ). Furthermore, Transwell assays were performed to examine how GCC2-AS1 affected the metastatic ability of LUAD cells. Depletion of GCC2-AS1 levels in the A549 and H292 cell lines notably reduced the numbers of invaded cells ( Figure 3C ). These data showed that GCC2-AS1 could enhance the proliferation and invasion of LUAD cells.

The “DESeq2” R package was used to perform DEGs analyses of LUAD samples from the high and low GCC2-AS1 expression groups as identified in the TCGA database. Based the threshold value (|log(FC)|>1.5 and adjusted P-value <0.05), a total of 361 genes (384 upregulated and 107 downregulated) were identified as being differentially expressed ( Figure 4A ). The top 10 genes with positive and negative correlations in the low and high GCC2-AS1 expression groups are shown in Figure 4B .

The R package clusterProfiler was used to performed a GO enrichment analysis of GCC2-AS1-related DEGs involved in several BPs, MFs, and CCs. The results ( Figure 4C , Supplementary Table 1 ) showed that the genes for BPs (81 terms) were mainly involved with digestion, cellular glucuronidation, and regulation of respiratory gaseous exchange. The genes for MFs (three terms) and CC (one term) were mainly involved with glucuronosyltransferase activity, extracellular matrix structural constituent, DNA-binding transcription activator activity, RNA polymerase II-specific, and Golgi lumen, respectively. Similarity, a KEGG (13 terms, Supplementary Table 2 ) enrichment analysis revealed that the DEGs were predominately related to chemical carcinogenesis, gluconeogenesis, and the peroxisome proliferators-activated receptors (PPAR) signaling pathway ( Figure 4D , Supplementary Table 2 ).

Next, the R package clusterProfiler was used to perform a GSEA hallmark analysis to identify genes with critical involvement in signaling pathways in the low and high GCC2-AS1 expression groups of LUAD patients. All the enriched key pathways were identified according to significant differences (|NES|>1, P-value <0.05, and FDR q value <0.25) in the enrichment of genes in the MSigDB collection (C2.all.v6.2.symbols), including lung cancer poor survival ( Figure 4E ), metastasis ( Figure 4F ), epithelial mesenchymal transition (EMT, Figure 4G ), proliferation ( Figure 4H ), and tumorigenesis ( Figure 4I ) in the high GCC2-AS1 expression phenotype, suggesting that GCC2-AS1 might play a role in LUAD progression.

The survmier R package was used to perform a Kaplan-Meier analysis that evaluated the prognostic value of GCC2-AS1 in the GCC2-AS1-high and GCC2-AS1-low groups. The results showed that patients with high levels of GCC2-AS1 expression had shorter OS time than patients in the GCC2-AS1-low groups ( Figure 5A ). Univariate and multivariate Cox regression analyses were performed to further investigate factors that correlated with a patient’s prognosis ( Table 3 ). The univariate analysis showed that T stage (HR = 1.668 [1.184–2.349], P = 0.003), N stage (HR = 2.606 [1.939–3.503], P<0.001), M stage (HR = 2.111 [1.232–3.616], P = 0.007), pathologic stage (HR = 2.975 [2.188–4.045], P<0.001), primary therapy outcome (HR = 2.818 [2.004–3.963], P<0.001), residual tumor (HR = 3.973 [2.217–7.120], P<0.001), tumor status (HR = 6.211 [4.258–9.059], P<0.001), and GCC2-AS1 expression (HR = 1.560 [1.164–2.092], P = 0.003) were significantly correlated with the OS time of LUAD patients ( Table 2 ). A multivariate analysis revealed that GCC2-AS1 expression (HR = 1.816 [1.113–2.964], P = 0.017), primary therapy outcome [(HR = 2.442 [1.401–4.254], P = 0.002), and tumor status (HR = 6.028 [3.339–10.884], P<0.001) independently predicted overall time. These data suggest that GCC2-AS1 could serve as an independent prognostic factor for LUAD patients.

Furthermore, a nomogram that integrated GCC2-AS1 with other independent prognostic factors identified by the multivariate Cox analysis was created to evaluate the probability of accurately predicting a prognosis. As shown in Figure 5B , variables that contributed to the total score for each case were used to determine the probability of a LUAD patient’s survival at 1, 3, and 5 years. The C-index for the nomogram was 0.758 (95% CI: 0.733–0.782). The calibration curve was very close to the actual curve (45-degree line, Figure 5C ), indicating strong agreement between the predicted probability and actual probability. In short, these findings indicated that the nomogram could provide a relatively accurate prediction of LUAD patient’s survival time.

To better understand the factors that affect the prognosis of LUAD patients, we conducted stratification analyses by creating survival charts and a forest plot. These analyses showed that a high level of GCC2-AS1 expression was strongly correlated with shorter OS time among patients with T stage ( Figure 6A ), N stage ( Figure 6B ), M stage ( Figure 6C ), pathological stage ( Figure 6D ), tumor status ( Figure 6E ), TP53 status ( Figure 6F ), and KRAS status ( Figure 6G ). In addition, the forest plot based on the cox regression analysis showed that GCC2-AS1 expression levels were of prognostic significance in the N0 subgroup (HR = 1.829 [1.190–2.813], P = 0.006) of N stage, stage I & stage II subgroups (HR = 1.496 [1.042–2.149], P = 0.029) of pathological stage, WT group (HR = 1.571 [1.111–2.220], P = 0.011) of KRAS status, Mut subgroup (HR = 1.716 [1.090–2.702], P = 0.020) of TP53 status, tumor subgroup (HR = 1.722 [1.179–2.514], P = 0.005) of tumor status, T2 subgroup (HR = 1.592 [1.077–2.354], P = 0.020) of T stage, and M0 subgroup (HR = 1.456 [1.029–2.060], P = 0.034) of M stage ( Figure 6H ).