Long non-coding RNA LINC02613 is a prognostic biomarker for breast cancer and correlates with the cell cycle and immune infiltration based on TCGA data.

Background: The development of transcriptomics has provided a new perspective to understand cancers. However, the role of long non-coding RNAs (lncRNAs) has not been fully elucidated. In the field of breast cancer, exploring the biological role and impact of lncRNAs on patient prognosis is of great importance.
Methods: Expression data of the lncRNA, LINC02613, were downloaded from The Cancer Genome Atlas (TCGA) dataset. The Kruskal-Wallis rank sum test, Wilcoxon rank sum test and logistic regression analysis were used to investigate the relationship between LINC02613 and clinicopathological characteristics. Cox regression and Kaplan-Meier analyses were used to assess the prognostic value of LINC02613. Gene ontology (GO) analysis, gene set enrichment analysis (GSEA) and immune infiltration analysis were used to explore the potential biological mechanisms of LINC02613.
Results: The expression of LINC02613 was downregulated in tumor samples and was related to poor overall survival (OS) in breast cancer patients. Statistical analysis revealed that low expression of LINC02613 was associated with Asian ethnicity, older age, higher histological grade and higher expression of estrogen receptor (ER), progesterone receptor (PR) and human epidermal growth factor receptor 2 (HER2). Positive expression of immunohistochemical markers (ER, PR and HER2), high pathological stage and N stage were independent predictors for low LINC02613 expression. Decreased expression of LINC02613 caused an increase in Wnt2 and Wnt3 protein expression, which activated the Wnt signaling pathway, thereby promoting cell cycle progression, DNA replication and glycolysis. At the same time, the antitumor immunity of the body was also weakened when the expression of LINC02613 was decreased. Conclusions: Low expression of LINC02613 predicts poor OS in breast cancer patients. LINC02613 is involved in the regulation of the cell cycle, DNA replication and glycolysis via the Wnt signaling pathway, and it is related to antitumor immunity. 2022 Translational Cancer Research. All rights reserved.

Introduction

Breast cancer is the most diagnosed malignant tumor and one of the leading causes of cancer-related mortality among females worldwide. In 2020, the number of newly diagnosed breast cancer cases was 2.3 million, representing 11.7% of all cancer cases. With 690,000 deaths caused by breast cancer, it is the fifth leading cause of cancer-related mortality (1). From 2017 to 2019, China was one of the top 3 countries in the world with the highest breast cancer incidence and deaths (2,3).

In clinical settings, immunohistochemical markers [estrogen receptor (ER), progesterone receptor (PR), human epidermal growth factor receptor 2 (HER2) and Ki-67] are used to classify patients into 4 different categories (luminal A, luminal B, triple negative breast cancer (TNBC) and Her2-enriched), which serve as a basis for decision making for surgery and adjuvant therapies. Although this strategy has improved patient prognosis to some extent in recent decades, it is still far from benefiting every patient (4). According to a study in 2019, the global mortality rate of breast cancer continued to show a slow upward trend from 1990 to 2015 with recurrence and metastasis remaining the leading cause of death among patients (5). Advancements in second-generation sequencing technology and bioinformatics technology have allowed continual exploration of the biological behavior of tumor cells and their relationship with immune cells in the tumor microenvironment (TME) as well as searches for additional potential regulatory targets to achieve more precise and individualized medical care.

Long intergenic non-coding RNAs (lncRNAs) are defined as RNA transcripts composed of more than 200 nucleotides. As nomenclature has evolved, they are simply labeled lncRNAs (6). The omitted term, ‘intergenic’, indicates their origin as follows: transcribed from regions of the genome that do not encode any proteins (7,8). LncRNAs were once considered useless because they lack open-reading frames and their expression is much lower than that of protein-encoding genes (9). However, the importance of lncRNAs is gradually being recognized as they play important roles in a series of biological processes, including gene expression control, scaffold formation and epigenetic control. LncRNAs can encode some short polypeptides (10,11). In addition, their dysregulation has been shown to be related to tumorigenesis (12).

LINC02613 is transcribed from the 2p22.1 region of the human genome. In this study, we downloaded the expression data of LINC02613 from The Cancer Genome Atlas (TCGA) dataset and explored its impact on the prognosis of breast cancer patients. On the basis of differentially expressed genes (DEGs), we demonstrated that LINC02613 affected the Wnt signaling pathway by influencing the expression of Wnt2 and Wnt3, which in turn regulated the cell cycle and DNA replication as well as promoted glycolysis, and we also demonstrated that LINC02613 was associated with antitumor immunity. We present the following article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-21-2479/rc).

Methods

Acquisition and preprocessing of RNA-sequencing (RNA-seq) data

We downloaded RNA-seq data in HTSeq-fragments per kilobase million (FPKM) format from UCSC XENA (https://xenabrowser.net/datapages/). Data without clinical information were omitted. Data in FPKM format were transformed then into TPM format to perform the analysis. Our study was in accordance with the publication guidelines provided by TCGA (https://cancergenome.nih.gov/publications/publicationguidelines) and did not contain any experiments involving human participants or animals performed by any of the authors. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

DEGs and functional enrichment analysis

All tumor samples were divided into two groups (high expression group and low expression group) according to the median expression level of LINC02613. The high expression group was labeled the test group, and the low expression group was labeled the control group. The DESeq2 package was used to compare expression data between the two groups (13). |Log fold change (logFC)| >2 and adjusted P<0.05 were regarded as the thresholds for DEGs. Functional enrichment analysis was performed using the Metascape database and its online analysis tools (http://metascape.org) (14).

Gene set enrichment analysis (GSEA)

GSEA is a computational method for verifying the extent to which a previously defined gene set is consistent with the expressed gene set. In the present study, signaling pathways with significant differences between the high and low LINC02613 expression groups were identified. The c2.all.v7.0.symbols.gmt curated from the Molecular Signatures Database (MSigDB) collection was used as the reference gene set. Permutations of gene sets were performed 1,000 times for each analysis. A false discovery rate (FDR) <0.25 and adjusted P<0.05 were regarded as statistically significant enrichments.

Immune infiltration analysis

Marker genes of immune cells, including dendritic cells (DCs), T cells, T follicular helper cells, cytotoxic cells, effector memory T cells, CD8+ T cells, Th1 cells, plasmacytoid DCs (pDCs), macrophages, activated DCs (aDCs), neutrophils, immature DCs (iDCs), NK CD56dim cells, NK cells, central memory T cells, gamma delta T cells, Tregs, mast cells, T helper cells, NK CD56bright cells, Th17 cells, eosinophils, Th2 cells and B cells, were obtained from a previously published paper (15). The single-sample Gene Set Enrichment Analysis (ssGSEA) method in the GSVA package was used to perform immune infiltration analysis by checking the expression level of these genes in published gene lists. Spearman correlation was then used to investigate the relationship between LINC02613 and immune cells.

Statistical analysis

The Wilcoxon rank sum test and Wilcoxon signed-rank test were used to compare LINC02613 expression between tumor samples and paratumor samples. The Kruskal-Wallis rank sum test and Wilcoxon rank sum test were used to analyze the relationship between multiple clinical features and the expression of LINC02613. Survival curves were constructed by the Kaplan-Meier method to assess the prognostic value of LINC02613 in breast cancer patients with overall survival (OS) (period from the date of randomization to the date of death or the last follow-up) as the primary endpoint as well as disease-specific survival (DSS) (period from the date of randomization to the date of death caused by breast cancer) and progression-free survival (PFS) (period from the date of randomization to the date of tumor progression or death from any cause) as the secondary end points. Cox regression models were used to determine independent prognostic factors. Variables with P<0.1 in univariate analysis were included in subsequent multivariate analysis. Hazard ratios (HRs) with 95% confidence intervals (CIs) were measured in all tests. All hypothesis tests were two-tailed, and a P<0.05 was considered statistically significant. Receiver operating characteristic (ROC) analysis was used to evaluate the discrimination ability of LINC02613 (16). All statistical analyses were performed in R (v.3.6.2).

Construction and assessment of nomogram

Based on the independent prognostic factors established by Cox regression analysis, a nomogram was constructed to show the predicted survival probabilities for 1-, 3- and 5-year OS using the rms package (https://cran.rproject.org/web/packages/rms/index.html). A calibration plot was then constructed to assess the predictive efficiency of the nomogram.

Results

Expression of LINC02613 is downregulated in tumors

After preliminary analysis, LINC02613 was found to be expressed at low levels in 24 tumor subtypes, including bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma (CESC), endocervical adenocarcinoma, colon adenocarcinoma (COAD), diffuse large B-cell lymphoma (DLBC), esophageal carcinoma (ESCA), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), acute myeloid leukemia (LAML), low-grade glioma (LGG) of the brain, liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), ovarian serous cystadenocarcinoma (OV), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), testicular germ cell tumors (TGCT), thyroid carcinoma (THCA), thymoma (THYM), uterine corpus endometrial carcinoma (UCEC) and uterine carcinosarcoma (UCS) (). We then explored the expression of LINC02613 in breast cancers. The Wilcoxon rank sum test and Wilcoxon signed rank test confirmed that LINC02613 expression was downregulated in breast cancer tissues compared to normal tissues and adjacent noncancerous tissues (). In addition, ROC analysis was performed to measure the discriminative power of LINC02613 in breast cancer. The area under the curve (AUC) was 0.925 (95% CI: 0.908–0.941), representing an efficient power of LINC012613 to discriminate tumors from normal tissues ().

Figure 1

Pancancer analysis of LINC02613. ns, P≥0.05; *, P<0.05; **, P<0.01; ***, P<0.001. TCGA, The Cancer Genome Atlas.

Figure 2

LINC02613 is downregulated in breast cancer tissues compared to normal and adjacent noncancerous tissues. (A) The expression of LINC02613 was significantly lower in breast cancer tissues than in normal tissues according to the Wilcoxon rank sum test (P<0.001); (B) LINC02613 showed significantly lower expression in breast cancer tissues than in adjacent noncancerous tissues according to the Wilcoxon signed-rank test (P<0.001).

Figure 3

Receiver operating characteristic (ROC) curve showing the discrimination power of LINC02613. The AUC was 0.925, suggesting an efficient power to discriminate tumors from normal tissues. AUC, area under the curve; TPR, true positive rate; FPR, false positive rate.

Relationship between LINC02613 and clinicopathological characteristics

Clinicopathological characteristics were compared between patients in the high expression and low expression groups (). A total of 1,065 patients were enrolled in our study, and 976 of whom had racial information, including 60 Asian patients, 179 African American patients and 737 Caucasian patients. Statistical results revealed that a greater proportion of Asian patients expressed lower LINC02613 than others (P=0.006). The mean age of the low expression group was 61 years (ranging from 50 to 70 years), while the mean age of the high expression group was 55 years (ranging from 47 to 64 years) with a statistically significant difference (P<0.001). Among all patients, those with left-sided neoplasms accounted for 51.9%, and those with right-sided neoplasms accounted for 48.1%. For pathological stage, lower expression of LINC02613 was related to higher pathological stage (stages III and IV vs. stages I and II, P=0.016). In addition, there were more patients with low expression of LINC02613 in the T4 stage and N2 stage (71.4% for T4, P=0.04; and 63.8% for N2, P=0.01). Low expression of LINC02613 was related to patients with positive expression of ER, PR and HER2 (ER positive vs. ER negative: 54.9% vs. 31.6%, P<0.001; PR positive vs. PR negative: 58.2% vs. 38.8%, P<0.001; HER2 positive vs. HER2 negative: 66.2% vs. 44.5%, P<0.001). The above results are shown in Figure S1. Further logistic analysis revealed that positive expression of ER [odds ratio (OR) =0.38; 95% CI: 0.28–0.52; P<0.001], PR (OR =0.52; 95% CI: 0.40–0.68; P<0.001) and HER2 (OR =0.41; 95% CI: 0.28–0.59; P<0.001), higher pathological stage (OR =0.71; 95% CI: 0.53–0.94; P=0.016) and N stage (OR =0.78; 95% CI: 0.61–1.00; P=0.047) were independent predictors for low expression of LINC02613 ().

Table 1

Characteristics of patients

Characteristics	Level	Low expression group (n=533)	High expression group (n=532)	P
T stage, n (%)	T1	130 (47.3)	145 (52.7)	0.04
	T2	313 (50.9)	302 (49.1)
	T3	63 (46.0)	74 (54.0)
	T4	25 (71.4)	10 (28.6)
N stage, n (%)	N0	236 (46.5)	271 (53.5)	0.01
	N1	175 (50.9)	174 (49.1)
	N2	74 (63.8)	42 (36.2)
	N3	35 (47.3)	39 (52.7)
M stage, n (%)	M0	449 (50.5)	440 (49.5)	1.000
	M1	10 (50.0)	10 (50.0)
Pathological stage, n (%)	Stage I & II	374 (47.6)	412 (52.4)	0.016
	Stage III & IV	144 (56.3)	112 (43.7)
PR status, n (%)	Negative	131 (38.8)	207 (61.2)	<0.001
	Positive	370 (58.2)	304 (41.8)
ER status, n (%)	Negative	75 (31.6)	162 (68.4)	<0.001
	Positive	427 (54.9)	351 (45.1)
HER2 status, n (%)	Negative	244 (44.5)	304 (55.5)	<0.001
	Positive	104 (66.2)	53 (43.8)
Race, n (%)	Asian	41 (68.3)	19 (31.7.)	0.006
	Black or African American	82 (17.4)	97 (19.2)
	White	349 (73.9)	388 (77.0)
Anatomic neoplasm subdivisions, n (%)	Left	275 (49.7)	278 (50.3)	0.877
	Right	258 (50.4)	254 (49.6)
Age, median (IQR), years		61.00 (50.00, 70.00)	55.50 (47.00, 64.00)	<0.001

P<0.05 indicates statistical significance. ER, estrogen receptor; PR, progesterone hormone receptor; HER2, human epidermal growth factor receptor 2; T, topography distribution; N, lymph node metastasis; M, distant metastasis; IQR, interquartile range.

Table 2

Results of logistic regression analysis

Characteristics	OR in LINC02613 expression	OR	P value
T stage (T3 & T4 vs. T1 & T2)	1,062	0.95 (0.68–1.31)	0.739
N stage (non-N0 vs. N0)	1,046	0.78 (0.61–1.00)	0.047
M stage (M1 vs. M0)	909	1.02 (0.41–2.51)	0.964
Pathological stage (stage III & stage IV vs. stage I & stage II)	1,042	0.71 (0.53–0.94)	0.016
PR status (positive vs. negative)	1,012	0.52 (0.40–0.68)	<0.001
ER status (positive vs. negative)	1,015	0.38 (0.28–0.52)	<0.001
HER2 status (positive vs. negative)	705	0.41 (0.28–0.59)	<0.001

P<0.05 indicates statistical significance. OR, odds ratio; ER, estrogen receptor; PR, progesterone hormone receptor; HER2, human epidermal growth factor receptor 2; T, topography distribution; N, lymph node metastasis; M, distant metastasis.

Low expression of LINC02613 indicates poor OS in breast cancer patients

Kaplan-Meier analysis revealed that low expression of LINC02613, although having no effect on DSS and PFS, predicted worse OS in patients [HR =0.65 (0.49–0.90), P=0.01] (). Subgroup analysis of the TNM staging system revealed that low LINC02613 expression was significantly associated with poor OS in patients with stage T2 disease [HR =0.581 (0.364–0.928), P=0.023], stage N0 disease [HR =0.540 (0.297–0.982), P=0.043] and stage M0 disease [HR =0.660 (0.459–0.949), P=0.025] (). T2 indicates that the maximum diameter of the primary tumor was >2.0 cm and <5.0 cm. N0 indicates that there was no infiltration of regional lymph nodes, and M0 indicates that no distant metastasis was discovered.

Figure 4

Kaplan-Meier analysis. (A) The Kaplan-Meier plot of relationship between LINC02613 and OS; (B) The Kaplan-Meier plot of relationship between LINC02613 and DSS; (C) The Kaplan-Meier plot of relationship between LINC02613 and PFS. OS, overall survival; DSS, disease-specific survival; PFS, progression-free survival.

Figure 5

Subgroup analysis revealing that low LINC02613 expression is significantly associated with poor OS of patients in the stage T2, stage N0 and stage M0 groups. If the HR of a group was infinite, the result was removed from the plot. OS, overall survival.

Prognostic factors identified by Cox regression analysis

Six independent prognostic factors were identified by univariate Cox regression, including T stage [HR =1.673 (1.152–2.429), P=0.007), N stage [HR =2.145 (1.497–3.073), P<0.001], M stage [HR =4.327 (2.508–7.465), P<0.001], pathological stage [HR =2.519 (1.787–3.549), P<0.001], age [HR =2.036 (1.468–2.822), P<0.001] and LINC02613 [HR =0.650 (0.469–0.901), P=0.010]. ER status (P=0.062) and HER2 status (P=0.059) were marginally significant. Further multivariate analysis confirmed that five of these factors, including T stage [HR =2.508 (1.293–4.864), P=0.007], M stage [HR =2.939 (1.141–7.545), P=0.026], ER status [HR =2.738 (1.597–4.694), P<0.001], age [HR =3.326 (2.001–5.501), P<0.001] and LINC02613 [HR =0.572 (0.332–0.984), P=0.544], were independent prognostic factors ().

Table 3

Results of Cox regression analysis

Characteristics	Total (N)	Univariate analysis HR (95% CI)	P value	Multivariate analysis HR (95% CI)	P value
T stage (T3 & T4 vs. T1 & T2)	1,061	1.673 (1.152–2.429)	0.007	2.508 (1.293–4.864)	0.007
N stage (N1 & N2 & N3 vs. N0)	1,045	2.145 (1.497–3.073)	<0.001	1.304 (0.699–2.435)	0.404
M stage (M1 vs. M0)	909	4.327 (2.508–7.465)	<0.001	2.934 (1.141–7.545)	0.026
Pathological stage (stage III & stage IV vs. stage I & stage II)	1,041	2.519 (1.787–3.549)	<0.001	1.953 (0.932–4.095)	0.076
PR status (negative vs. positive)	1,011	1.312 (0.931–1.849)	0.120
ER status (negative vs. positive)	1,014	1.420 (0.983–2.052)	0.062	2.738 (1.597–4.694)	<0.001
HER2 status (negative vs. positive)	705	0.621 (0.378–1.019)	0.059	1.210 (0.657–2.227)	0.541
Age (>60 vs. ≤60 years)	1,064	2.036 (1.468–2.822)	<0.001	3.326 (2.011–5.501)	<0.001
Race (Caucasian vs. Asian and African American)	975	0.880 (0.593–1.306)	0.526
Anatomic neoplasm subdivisions (right vs. left)	1,064	0.776 (0.559–1.077)	0.130
LINC02613 (high vs. low)	1,064	0.650 (0.469–0.901)	0.010	0.572 (0.332–0.984)	0.044

P<0.05 indicates statistical significance. Variables with P<0.1 in univariate analysis were included in multivariate analysis. ER, estrogen receptor; PR, progesterone hormone receptor; HER2, human epidermal growth factor receptor 2; T, topography distribution; N, lymph node metastasis; M, distant metastasis; HR, Hazard ratio.

Construction and evaluation of the nomogram

A nomogram was constructed to provide surgeons with a quantitative approach to predict the prognosis of breast cancer patients (). In this nomogram, each variable corresponds to a scaled line segment, in which length represents the rescaled range and impact on prognosis. The individual score of each variable and their sum can be easily obtained, and the probability of survival at different time nodes can be obtained. A calibration plot was then constructed to intuitively evaluate the prediction performance of the nomogram (). The abscissa represents the probability of prognosis obtained from the nomogram, and the ordinate represents the actual prognosis (0 to 1 indicates that the probability of events is 0% to 100%). The colored lines are fitted lines, and the diagonal represents the ideal situation. Increased overlapping indicated better prediction performance. The C-index of the nomogram was 0.713 (0.687–0.738), indicating a medium level of predictive accuracy.

Figure 6

Nomogram for predicting 1-, 3-, and 5-year OS for breast cancer patients. The top row shows the point value for each variable. Rows 2–6 indicate each variable included in this nomogram. Row 7 shows total points of all 5 variables. Row 9–11 shows corresponding survival probability.

Figure 7

Calibration plot of the nomogram for predicting 1-, 3-, and 5-year OS for breast cancer patients. OS, overall survival.

Identification of DEGs and gene ontology (GO) analysis

Because the importance of LINC02613 for breast cancer patients has been well established, we were interested in further exploring its functions in organisms. Differential expression analysis revealed that the expression of 382 genes met the threshold of DEGs, of which 126 were upregulated and 256 were downregulated. Ten of these genes were selected for heatmapping (). GO analysis revealed that genes regulated by LINC02613 were mainly enriched in biological processes related to tissue development and energy metabolism, including cornification, regulation of neuroblast proliferation, pentose and glucuronate interconversions and galactose metabolism. The detailed results and a network of their interactions are provided in Figure S2.

Figure 8

Identification of DEGs. (A) Volcano plot of LINC02613-related DEGs; (B) heatmap of LINC02613-related DEGs. DEGs, differentially expressed genes.

GSEA

Traditional functional enrichment analysis only focuses on a few genes with significant differential expression, omitting the genes with great biological importance but inconspicuously different expression. Importantly, GSEA fills this gap. The GSEA results showed that six datasets were significantly enriched in the LINC02613 low expression group, including DNA replication, cell cycle, glycolysis gluconeogenesis, antigen processing and presentation, chemokine signaling pathway and WNT signaling pathway (Figure S3). Tumor cells are characterized by their unrestricted ability to divide and proliferate; DNA replication and smooth cell cycle are necessary processes for cell proliferation. In addition, adequate nutrition is a prerequisite, but the energy provided by normal intracellular aerobic oxidation is insufficient to support this abnormal process. According to Warburg's theory, in such a condition, tumor cells will generate energy by activating glycolysis in the presence of oxygen (17). Antigen processing and presentation as well as chemokine secretion are important components of antitumor immunity (18). Combined with previous studies, we hypothesized that low expression of LINC02613 regulates tumor cell cycle progression, DNA replication and glycolysis through the Wnt signaling pathway, and we also hypothesized that low expression of LINC02613 weakens antitumor immunity.

Relationship between LINC02613 and cell cycle

The cell cycle refers to the activity of a cell from the beginning of one mitosis to the end of the next mitosis. The cell cycle is divided into the following five phases: G1 phase (pre-DNA synthesis), S phase (DNA synthesis), G2 phase (post-DNA synthesis), M phase (mitosis) and G0 phase (pre-DNA synthesis). DNA replication occurs in S phase. When growth factors bind to cell surface receptors, cells in the G0 phase are activated and express cyclins and cyclin-dependent kinases (CDKs). Cyclin D1, the first discovered and most studied cyclin molecule, is a cell cycle initiator and growth-sensing factor. In the early G1 phase, cyclin D1 binds to CDK4 or CDK6 to drive cell cycle progression and DNA replication (19). In this study, we examined the relationship of CCND1 (the gene encoding cyclin D1), CDK4 and CDK6 with LINC02613. The results showed that the expression of CCND1 was negatively correlated with LINC02613 (Figure S4).

Relationship between LINC02613 and the Wnt signaling pathway

The Wnt pathway is an evolutionarily conserved signaling pathway, and it can be divided into β-catenin-dependent and β-catenin-independent subtypes. Aberrant activation of the β-catenin-dependent Wnt pathway has been shown to result in tumorigenesis and increased cellular glycolysis (20,21). When Wnt proteins bind to Frizzled receptors, repressed β-catenin is released, activating downstream TEF/LEF and affecting the transcription of numerous genes, including CCND1. There are many members of the Wnt protein family. Wnt2, Wnt3 and Wnt3a are all considered to be initiators of the β-catenin-dependent Wnt signaling pathway. In this study, we investigated the relationship between the Wnt2, Wnt3 and Wnt3a coding genes and LINC02613. The results revealed that the expression of Wnt2 and Wnt3 was negatively correlated with LINC02613 (Figure S5). Considering the relationship between LINC02613 and CCND1, the downregulation of LINC02613 may promote the expression of Wnt2 and Wnt3, leading to activation of the Wnt signaling pathway, which promotes the expression of downstream cyclin D1, leading to cell cycle progression and DNA replication of tumor cells.

Relationship between LINC02613 and immune infiltration

The balance between the tumoricidal effect of immune cells and the immune escape of tumor cells in the microenvironment plays an important role in the occurrence and metastasis of cancer. We analyzed the relationship between the expression of LINC02613 and infiltration of immune cells in the TME by ssGSEA. The results suggested that the expression of LINC02613 was positively related to the infiltration of 20 out of 24 immune cells, including CD8+ T cells, NK cells, macrophages and DCs, which remain the workhorses in eliminating tumors. In addition, other innate immune cells, such as mast cells and neutrophils, as well as functional T cells, including Tfh, Tem, Th1, Tcm, Tgd, Treg and T helper cells, were also shown to be highly expressed in the test group (). These results suggested that the expression of LINC02613 is positively related to immune infiltration in the TME and that the downregulation of LINC02613 may weaken the antitumor effect in breast cancer patients.

Figure 9

Relationship between LINC02613 and immune infiltration. (A) Immune infiltration in test group; (B) LINC02613 was positively related to CD8+ T cells; (C) LINC02613 was positively related to macrophages; (D) LINC02613 was positively related to NK cells; (E) LINC02613 was positively related to DCs. DCs, dendritic cells; pDCs, plasmacytoid DCs; aDCs, activated DCs; iDCs, immature DCs

Discussion

Treatment for breast cancer is entering the era of individual-based precision medicine. LncRNAs have great potential to stratify the risk of recurrence and predict the prognosis of patients (22).

Dysregulation of lncRNA expression has been shown to be related to the occurrence and metastasis of tumors. HOTAIR has been shown to be highly expressed in breast cancer tissues and to promote tumor cell metastasis by reprogramming the chromatin state (23). LUNAR has been shown to regulate T-acute lymphoblastic leukemia (T-ALL) through the Notch signaling pathway (24). Another non-coding RNA, MALAT1, has been shown to be associated with tumor cell metastasis in a variety of cancer types, including lung, breast and colon cancers (25). In the present study, we identified a potential breast cancer suppressor gene, LINC02613.

Analysis of TCGA database revealed that LINC02613 expression was significantly lower in tumor samples compared to normal samples. Low expression of LINC02613 was related to ethnicity (Asian vs. non-Asian, P<0.001), mean age (61.0 vs. 55.5 years, P<0.001), pathological stage (stages I and IV vs. stages I and II, P=0.012), T and N stage (T4 vs. non-T4, P=0.012; non-N0 vs. N0, P=0.038) and expression of immunohistochemical markers (ER, PR and HER2). Logistic regression analysis revealed that positive expression of immunohistochemical markers (ER, PR and HER2), higher pathological stage (stages III and IV vs. stages I and II, P=0.016) and N stage (non-N0 vs. N0, P=0.047) were independent predictors for low expression of LINC02613. These results indicated that LINC02613 may play a tumor suppressive role and that low expression of LINC02613 implies a more aggressive and metastatic ability of breast cancer cells.

Based on TCGA database, Kaplan-Meier survival analysis and Cox regression analysis indicated that low expression of LINC02613 was associated with worse OS in breast cancer patients. We identified the following five independent prognostic factors: age, T stage, LINC02613, metastasis and ER status. These results indicated that LINC02613 may serve as a candidate marker to assist survival prediction in clinical settings. Therefore, a nomogram including these prognostic factors was subsequently built and showed good predictive value.

Because the importance of LINC02613 for breast cancer patients has been well established, we further explored its functional mechanisms. In conventional GO analysis, LINC02613 was mainly associated with energy metabolism and tissue development. In GSEA, six gene sets were significantly enriched in the LINC02613 low expression group, including DNA replication, cell cycle, glycolysis gluconeogenesis, Wnt signaling pathway, antigen presentation and chemokine signaling pathway. Further gene-to-gene correlation analysis showed that LINC02613 expression was inversely correlated with Wnt2, Wnt3 and CCND1.

The eukaryotic cell cycle is a highly organized and strictly regulated process that ensures the correct distribution of genetic material in the nucleus during cell division (26). Cyclin D1, as one of the key regulators of the cell cycle, binds to CDK4 or CDK6, enters the nucleus and phosphorylates downstream molecules to promote the transcription of genes required for cell proliferation (27). CCND1, as the encoding gene of cyclin D1, has been shown to be overexpressed in tumor cells and is regarded as a negative prognostic biomarker (28). Increasing evidence has revealed that CCND1 is regulated by the Wnt signaling pathway (29,30). When the Wnt pathway is activated, the downstream TEF/LEF complex promotes the transcription of CCND1, which induces cell cycle progression and DNA replication (31). In the present study, the expression of Wnt2, Wnt3 and CCND1 was negatively correlated with LINC02613, which implied that low expression of LINC02613 contributed to the unrestricted proliferation of tumor cells by upregulating Wnt protein, thus activating the Wnt pathway and promoting transcription of downstream CCND1.

Metabolic reprogramming is a prominent feature of tumors, and abnormalities in glucose metabolism were the first identified examples of metabolic reprogramming (32,33). Under aerobic conditions, tumor cells still rely on glycolysis for energy, which is also known as the Warburg effect (17). Recent studies have revealed that the canonical Wnt signaling pathway promotes tumor cell glycolysis in a variety of ways, including activating the expression of MCT-1, CYC1 and ATP synthase through the downstream transcription factor, TCF/LEF, which promotes the overexpression of c-Myc and regulates the activity of rate-limiting enzymes in the energy metabolism pathway (21). In the present study, both glycolysis gluconeogenesis and the Wnt signaling pathway were enriched in the low LINC02613 expression group, implying that low expression of LINC02613 may promote cellular glycolysis by activating the Wnt pathway.

In addition, the present study showed that LINC02613 expression was positively correlated with several types of antitumor immune cells, including CD8+ T cells, macrophages, DCs and Treg cells, which are recognized to play vital roles in antitumor immune responses (34). It is conceivable that the number of immune cells in the TME decreases dramatically when LINC02613 is expressed at low levels, further exacerbating tumor growth and metastasis.

In the present study, we first explored the relationship between LINC02613 and breast cancer. Despite the important findings in the present study, there were several limitations. First, although advances in sequencing and database technologies have greatly broadened information access compared to the past, the sample size in this study was still limited due to the single source of sample information. Second, due to objective constraints, we only explored the expression of LINC02613 on the bioinformatics platform, which is insufficient for rigorous scientific study. Thus, there is a clear need for additional multicenter and large-scale clinical trials before applying our nomogram in clinical practice. In addition, further cytological and in vivo experiments are required to validate our findings in terms of molecular mechanisms.

Conclusions

LINC02613 is a potential tumor suppressor gene. Low expression of LINC02613 predicts poor OS in breast cancer patients. LINC02613 not only affects the Wnt signaling pathway by influencing the expression of Wnt2 and Wnt3 proteins, which in turn regulates the cell cycle, regulates DNA replication and promotes glycolysis, but is also associated with antitumor immunity.