First, the “TCGA-ALL-FPKM” dataset (sample=11093) was used to compare the differential expression of SLC6A8 in tumor tissues and in normal ones in 33 cancers using the Xiantao Academic Online website (https://www.xiantao.love/) (P-value<0.05; wilcoxon rank sum test). Moreover, NSCLC was selected as the study subject, and the transcriptional data of SLC6A8 in NSCLC (tumor=1037; normal=108) was downloaded from the TCGA database (https://portal.gdc.cancer.gov/). Differential expression of SLC6A8 in tumor and normal tissues was analyzed using the “limma” R package (P-value<0.05; t-test). Next, the cBioPortal online website (http://www.cbioportal.org/) was used to select 3 LAUD TCGA datasets (TCGA Firehose Legacy; TCGA PanCancer Altas; TCGA Nature 2014) and 3 LUSC TCGA datasets (TCGA Firehose Legacy; TCGA PanCancer Atlas; TCGA Nature 2012) to explore the genetic alterations of SLC6A8 in NSCLC. Finally, the “polts” section was selected to visualize the relationship between SLC6A8 expression and a copy number.

Pre-experimental tissue microarrays (TMA) 1 included 57 pairs of NSCLC tissues and paracancerous tissues (18 females and 39 males), purchased from Superbiotek (Shanghai, China). The mean age of patients participating in the study was 61.1 years (range: 34-84 years). (Stage: T1aN0M0 to T4N0M1c; 2004 World Health Organization criteria). TMA1 was mainly used to study paired difference analysis of SLC6A8 in NSCLC and in paraneoplastic tissues (P-value<0.05; t-test). TMA2 obtained 140 NSCLC tissue samples and 10 normal adjacent tissue samples from patients who underwent surgical resection in the Department of Thoracic Surgery, Zhongshan Hospital, Fudan University, from January 2005 to December 2005. All patients had complete clinical information (38 females and 112 males) and a mean age of 60.1 years (range, 26-79 years) for NSCLC (stage: Ia to IIIa; American Joint Committee on Cancer and the Union for International Cancer Control criteria). Follow-up of this clinical information was recorded until July 2013. TMA2 was used to analyze the unpaired difference expression of SLC6A8 in NSCLC tissues and normal paraneoplastic tissues (P-value<0.05; t-test). Neither TMA1 nor TMA2 patients received chemotherapy, radiotherapy or biological therapy before surgery.

To detect the expression of SLC6A8 in NSCLC, immunohistochemistry was performed using the standard indirect immunoperoxidase procedures. Paraffin specimens were cut 4-µm thick, mounted on slides, baked, deparaffinized and hydrated according to conventional methods. 200 ml of 3% H2O2 and 1 ml of NaN3 were used to inactivate the endogenous peroxidase activity, followed by antigen recovery, which was performed with 10 mM sodium citrate buffer (pH 6.0). Slides were incubated for 1 hour at room temperature in 10 mM TBS with 4% normal rabbit serum (Proteintech) and incubated with primary antibody against SLC6A8 (1:50, 20299-1-AP, Proteintech Group, Inc, China) at 4°C overnight. Thereafter, the secondary antibody (1:200, K5007, DAKO, China) was developed for 35 minutes at 37°C. Finally, the specimens were weakly re-stained with hematoxylin at 37°C, dehydrated and covered with coverslips.

To quantify the expression of SLC6A8 protein in NSCLC tissues, the image was converted to grayscale by selecting “8 bit” and the grayscale values were converted to optical density values using “Uncalibrated OD” on ImageJ (22) software. Using the “Set Measurement” module, the area to be stained on the image can be set. The formula for calculating the average optical density (AOD) is as follows:

Paired difference analysis of AOD in tumor tissues and normal tissues in TMA1 was performed using GraphPad Prism software, while unpaired difference analysis was performed in TMA2 (P-value<0.05; t-test). Finally, TMA2 was divided into high- and low-expression groups according to the median AOD value of SLC6A8, and the significance of survival between the two groups was analyzed using GraphPad Prism software (P-value<0.05; log-rank test).

To assess the correlation between SLC6A8 expression and clinicopathological characteristics, we analysed the correlation between SLC6A8 expression and T stage, N stage, M stage, pathological stage, gender, age and smoking using TCGA data selected from the Xiantao Academic Online website (P-value<0.05). The impact of these clinicopathological factors on overall survival was then analysed by means of univariate and multivariate COX analysis (P-value<0.05).

Firstly, the miRNA data for NSCLC was downloaded from the TCGA database by selecting the “miRNA Expression Qualification” section (cancer=999; normal=91). The miRNA expression matrix was obtained by organizing the raw data using perl software. The Starbase online database (http://starbase.sysu.edu.cn/), a rich database comprised of miRNA-ncRNA, miRNA-mRNA, RBP-RNA and RNA-RNA data that provides a variety of visual interfaces for exploring microRNA targets by searching for microRNA targets through high-throughput CLIP-Seq experimental data and degradome experimental data (23), was then used to find SLC6A8-binding miRNAs by selecting “miRNA-mRNA” in the “miRNA-Target” module (programNum>=2). Moreover, the obtained miRNA-mRNA data were subjected to co-expression network construction by Cytoscape software. Since miRNAs and mRNAs are negatively correlated in ceRNA networks, cor<-0.2 and P-value<0.001 were set as correlation filter conditions in order to find SLC6A8-related miRNAs in NSCLC. LogFC<0 and diffP-value<0.05 were filtered conditions for differential miRNA expression in tumor tissues and normal tissues.

To explore the targeting relationship between “miRNA-lncRNA”, the screened miRNAs were entered into the “miRNA-Target” module in Starbase, and the screened data were downloaded (programNum>=2). Then, the obtained miRNA-lncRNA data were used for co-expression network construction using Cytoscape software. Since miRNAs and lncRNAs were negatively correlated in the ceRNA network, cor<-0.2 and P-value<0.001 were used as correlation filtering conditions. LogFC>0 and diffP-value<0.05 were used as filtering conditions for differential lncRNA expression in tumor tissues and normal tissues. Finally, the ceRNA networks were obtained by visualizing the ceRNA networks using Cytoscape software.

Firstly, lncRNA data downloaded in TCGA were used to analyze the correlation of SLC6A8 with lncRNAs in the regulatory network (cor>0.2, P-value<0.001), and these data were used to further compare the differential expression of lncRNA in NSCLC tissues versus normal tissues (LogFC>0, diffP-value<0.05; t-test). According to the median values of miRNAs and lncRNAs expression in NSCLC, the miRNA and lncRNA samples downloaded from TCGA were divided into high and low expression groups, respectively, and the survival differences between the two groups were analyzed using the log-rank test (P-value< 0.05). Then, the sequences of each lncRNA were queried using the LNCipedia database (https://lncipedia.org/), and the obtained sequences were entered into the lncLocator database (http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/) to obtain the percentage distribution of each lncRNA in the cell. Finally, base pairing between prognosis-related lncRNA-miRNA and mRNA-miRNA were predicted using the Starbase online database.

First, SL6A8 methylation levels between LUAD and LUSC tissues and normal tissues were analyzed using the ENSG00000130821 dataset selected from the EWAS Data Hub online database (https://ngdc.cncb.ac.cn/ewas/datahub/index). Thereafter, the association between the expression of SLC6A8 and their DNA methylation status was investigated using MEXPRESS (https://mexpress.be). Finally, the prognosis of SLC6A8-associated methylation sites in NSCLC was analyzed through the EWAS Data Hub online database.

The correlation between SLC6A8 expression and marker gene expression on the surface of six immune cells in NSCLC, the difference between 79 immunomodulators in the high and low SLC6A8 expression groups and the correlation between 79 immunomodulators and SLC6A8 expression were first analyzed by the “limma” R package using TCGA data (P-value<0.05). Then, the “survival” R package was used to screen for prognosis-related immunomodulator genes (P-value<0.05). To explore the immuno-prognostic characteristics of SLC6A8 in NSCLC, we selected prognosis-related immunomodulator genes (cor>0.2) to construct a COX risk proportional regression model. The formula for calculating the risk score was as follows:

The “coefi” and “Xi” represent the coefficient and expression level of each prognosis-related immunomodulator, respectively. Classification of NSCLC patients into high- and low-risk groups was based on median model scores, and a comparison of survival differences between the two groups using the “survivor” and “survminer” R packages (log-rank test; P-value<0.05) was made. Then, a combination of model risk scoring with clinical factors (age, gender, and stage) was made to explore whether risk scoring is an independent prognostic factor for NSCLC. Finally, the efficacy of immunotherapy was assessed between high- and low-risk groups through the Tumor Immune Dysfunction and Exclusion (TIDE) online website. (http://tide.dfci.harvard.edu/).

This study was reviewed and approved by the Research Ethics Committee of The Affiliated Changzhou No. 2 People’s Hospital of Nanjing Medical University. All patients were informed and voluntarily signed an informed consent form for the collection of clinical tissue samples. All specimens were processed and anonymized according to ethical and legal standards.

We implemented all statistical analyses with R (version 4.1.0) and GraphPad Prism 9 (version 9.1.2). The difference in SLC6A8 expression between the two groups was compared by t test, and the difference in survival between the two groups was compared by log-rank test. P-value<0.05 was considered statistically significant.

The results revealed that SLC6A8 was differentially expressed in various cancers, with SLC6A8 expression significantly higher in LUAD and LUSC than in normal tissues ( Figure 1A ; P-value<0.05). Analysis of the data downloaded from TCGA showed that SLC6A8 was significantly overexpressed in NSCLC than in normal tissues ( Figure 1B ; P-value<0.05) or in paraneoplastic ones ( Figure 1C ; P-value<0.05). Genetic alterations analysis results indicated three forms of alterations in SLC6A8: amplification, mutation and deep deletion, with gene amplification being the most common, followed by gene mutation ( Figure 1D ). In addition, the results of the copy number analysis also demonstrated that NSCLC samples with the SLC6A8 amplification exhibited higher mRNA expression ( Figure 1E ).

To prove the above findings, immunohistochemical staining experiments were performed to analyze the paired difference analysis of SLC6A8 between 57 pairs of NSCLC and paracancerous tissues and the unpaired difference analysis between 140 tumor samples and 10 normal samples. The results demonstrated that both paired difference and unpaired difference analysis of SLC6A8 in NSCLC showed statistical significance, and SLC6A8 expression was significantly higher in tumor tissues than in paraneoplastic ( Figure 2A–E ; P-value<0.05) or normal tissues ( Figure 2F–I ; P-value<0.05). Finally, TMA2 was integrated with data from clinical follow-up and divided into high- and low-expression groups based on the median value of the AOD of SLC6A8 (0.121), which revealed that the high-expression group showed a poorer prognosis compared to the low-expression group ( Figure 2J ; P-value<0.05).

Analysis of clinicopathological characteristics revealed significant differential expression of SLC6A8 in T stage, N stage, M stage, pathological stage, gender, age and smoking ( Table 1 ; P-value<0.05). In addition, univariate COX analysis showed that T stage (T1&T2 vs. T3&T4), N stage (N0&N1 vs. N2&N3), M stage (M0 vs. M1), pathological stage (Stage I vs. Stage II&III&IV) and age (<=65 vs. >65) were prognostically relevant risk factors. Multivariate COX analysis indicated that these factors could be independent risk factors for prognosis ( Table 2 ; P-value<0.05).

To explore the lncRNA-miRNA-mRNA networks that regulate the expression of SLC6A8 in NSCLC, the expression matrix of miRNAs from the TCGA database combined with SLC6A8-related miRNAs from the Starbase database (miRNA=41) was extracted to construct a miRNA-mRNA co-expression network ( Figure 3A ). In the miRNA-mRNA correlation analysis, there was a correlation between hsa-miR-26a-5p expression and SLC6A8 expression in NSCLC (cor<-0.2, P-value<0.001; Figure 3B ), and hsa-miR-26a-5p expression was significantly lower in NSCLC than in normal tissues (P-value<0.05; Figure 3C ). Analysis using the Starbase database identified 84 target lncRNAs (lncRNA=97) associated with hsa-miR-26a-5p ( Figure 3D ). Among them, LINC01703, AC104088.1, DLX6-AS1, AC013652.1 and AL513318.2 met the filtering conditions for correlation analysis with hsa-miR-26a-5p (cor<-0.2, P-value<0.001), and they were significantly overexpressed in NSCLC than in normal tissues (logFC<0, diffP-value<0.05) ( Table 3 ). Finally, the ceRNA networks that regulate SLC6A8 expression in NSCLC were mapped using Cytoscape software ( Figure 3E ).

Survival analysis revealed that the low expression group of hsa-miR-26a-5p in NSCLC had a poorer prognosis compared to the high expression group (cut point=11.68; Figure 4A ). Meanwhile, hsa-miR-26a-5p was significantly less expressed in NSCLC tissues than in normal tissues ( Figure 3C ), so hsa-miR-26a-5p could be used as a prognostic biomarker for NSCLC. The study demonstrated that LINC01703, AC104088.1, DLX6-AS1, AC013652.1 and AL513318.2 were significantly associated with SLC6A8 (cor>0.2, P-value<0.001) and were significantly differentially expressed in NSCLC tissues versus normal tissues (logFC<0, P-value< 0.05) ( Table 3 ). lncRNAs survival analysis showed that the AC104088.1 (cut point=1.172) and DLX6-AS1 (cut point=0.4318) low expression groups had a poor prognosis compared to the high expression group, while the AL513318.2 (cut point=0.5129) high expression group was associated with a poor prognosis ( Figure 4B ). Cellular localization analysis showed that AC104088.1, DLX6-AS1 and AL513318.2 had the highest percentage in the cytoplasm ( Figure 4C ). Therefore, AC104088.1, DLX6-AS1 and AL513318.2 can be used as prognostic biomarkers for NSCLC ( Figure 4D ). Since SLC6A8 is highly expressed in NSCLC and associated with a poor prognosis, the AL513318.2/hsa-miR-26a-5p/SLC6A8 regulatory network could be used as a biomarker for poor prognosis and a new target for the treatment of NSCLC ( Figure 4E ). In addition, base pairing between AL513318.2/hsa-miR-26a-5p and SLC6A8/hsa-miR-26a-5p were predicted using the Starbase online database ( Figure 4F ).

To further explain the mechanism of aberrant upregulation of SLC6A8 in NSCLC tissues, the correlation between SLC6A8 expression levels and its methylation status was explored. Analysis of the ENSG00000130821 dataset in the EWAS Data Hub online database showed that methylation levels in LUAD ( Figure 5A ; Supplementary Material 1 ) and LUSC ( Figure 5B ; Supplementary Material 2 ) negatively correlated with SLC6A8 expression, and the methylation levels in tumor tissues were lower than those in normal tissues. Then, the MEXPRESS database exhibited multiple methylation regulatory sites in LUAD ( Figure 5C ) and LUSC ( Figure 5D ) that were negatively associated with SLC6A8 expression (r<-0.2, P-value<0.05). Among them, cg18010669, cg16632275, cg04676446, cg18476631, cg04307491, cg15731001, cg01577514, cg20095274, cg20765985, cg12177562 and cg13762255 are involved in SLC6A8 expression regulation in both LUAD and LUSC. Methylation site survival analysis indicated that 11 methylation sites expression were significantly associated with LUAD and LUSC survival ( Supplementary Material 3 ). In addition, cg12177562, cg13762255 and cg20765985 methylation sites had significantly poorer OS in the low expression group compared to the high expression group in both lung adenocarcinoma and lung squamous carcinoma.As these sites were significantly negatively correlated with SLC6A8 expression in lung adenocarcinoma and lung squamous carcinoma, these three sites may be potential therapeutic targets for modulating the poor prognosis of SLC6A8 in NSCLC.

Immunological correlation analysis showed that SLC6A8 expression mainly correlated negatively with B cell, CD8+ T cell, M1 macrophage, M2 macrophage, Neutrophil and Dendritic cell surface marker genes, with the most negative correlation with NRP1 ( Table 4 ). Research indicated that immunomodulators are an important component of the tumor immune microenvironment, and their alterations are closely related to the prognosis of patients (24, 25). In this study, 54 of the 79 immunomodulators had significant differences between the high and low SLC6A8 expression groups ( Figure 6A ), of which 31 immunomodulators predominantly correlated negatively with SLC6A8 expression, with CD40LG being the most negatively correlated ( Table 5 ). Then, prognostic analysis of 31 SLC6A8-associated immunomodulators in NSCLC indicated BLTA, CD160, CD40LG, and TNFRSF13C as low-risk factors and NT5E as high-risk factors ( Figure 6B ). To explore SLC6A8-mediated immune prognosis, risk proportional regression models of SLC6A8 prognosis-associated immunomodulator genes ( Figures 6C, D ) were constructed. The model risk scoring formula is as follows: risk core = (coefficient × Expression CD40LG) + (coefficient × Expression NT5E) The high-risk group classified according to the median value (1.024) of the model risk scoring had poorer survival compared to the low-risk group ( Figure 6E ). Moreover, model risk scoring combined with clinical factors in univariate and multifactorial COX analyses suggested that age, stage and risk scoring could be independent risk factors for NSCLC prognosis ( Figure 6F ). Finally, the effect of immunotherapy was evaluated between the high- and low-risk groups using the TIDE database, and the results revealed that immunotherapy was significantly better in the high-risk group than in the low-risk group ( Figure 6G ).