TNF-like ligand 1A is associated with progression and prognosis of human gastric cancer.

PURPOSE: This study aimed to investigate the function of TNF-like ligand 1A (TL1A) in the tumorigenesis and progression of gastric cancer (GC).
METHODS: RNA-seq gene expression and clinical information for GC patients were obtained from The Cancer Genome Atlas (TCGA) database. Differentially expressed genes (DEGs) between GC tissue samples and normal controls were screened with the edgeR package. Identification of gene co-expression and functional enrichment analyses were performed with Pearson's correlation analysis and gene set enrichment analysis (GSEA), respectively. Lastly, survival analysis was performed using the Kaplan-Meier method with the log rank test.
RESULTS: TL1A expression in GC tissue samples were significantly higher than that in normal controls (LogFC=1.07 and P=8.90E-07). Moreover, 215 genes, co-expressed with TL1A, and 21 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were obtained. Next, the miRNA-lncRNA/mRNA network, comprising 7 miRNAs, 27 lncRNAs, and 21 mRNAs, was constructed based on key genes from intersections between co-expression analysis and GSEA. In addition, survival analysis results demonstrated that TL1A (P=2.6e-07) was significantly associated with the overall survival (OS) of GC patients.
CONCLUSION: TL1A was involved in the tumorigenesis and progression of GC, and was significantly associated with the OS of GC patients.

Introduction

Gastric cancer (GC) is the fifth most common malignancy, and the third leading cause of cancer-related deaths (723,100 in 2012) globally.1 Currently, a series of risk factors, such as dietary factors, Epstein-Barr virus, and Helicobacter pylori infection have been identified for GC. Moreover, an increasing amount of evidence show that genetic factors and their downstream effects on cellular processes play important roles in GC progression, but, the pathogenesis of GC still remains complex.2,3 Additionally, despite the advances in cancer diagnosis and treatment, the prognosis of GC patients still remains unsatisfactory with the 5-year survival rate less than 30% due to its aggressive behavior.4 Hence, it is urgent to identify novel biomarkers and effective therapeutic targets for patients with GC.

TNF-like ligand 1A (TL1A; also known as TNFSF15 or VEGI), is a unique cytokine produced by endothelial cells, T cells, macrophages, monocytes and dendritic cells, which inhibits proliferation and angiogenesis of endothelial cells.5,6 TL1A has been reported to be involved in modulating vascular homeostasis and inflammation.7 Previous studies have reported that TL1A plays an important role in a variety of autoimmune inflammatory diseases, such as rheumatic arthritis (RA), asthma and inflammatory bowel disease.8 Additionally, TL1A is involved in the development of multiple tumor types,9,10 and is recognized as the potential therapeutic target for cancer therapy.11,12 Recently, reduced TL1A expression was found to be associated with poor prognosis in breast cancer patients.11 However, the role of TL1A in GC progression remains unclear.

In the present study, in order to reveal the potential role of TL1A in tumorigenesis of GC, we screened the differentially expressed genes (DEGs) and co-expressed genes associated with TL1A in GC. Subsequently, we performed gene set enrichment analysis (GSEA) to obtain the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. Finally, we constructed a TL1A-related miRNA-lncRNA/mRNA network and predicted the prognostic value of TL1A in GC.

Materials and methods

Data collection

Clinical and survival data from 418 patients with GC were retrieved from The Cancer Genome Atlas (TCGA) through the Genomic Data Commons (GDC) portal. Further, 375 tumor samples and 50 normal tissues samples with the RNA-seq gene expression data (HTSeq–Counts) were included in the present study. All data in our study were downloaded from the database of University of California Santa Cruz (UCSC) xena (https://xenabrowser.net).

Data preprocessing and differential expression analysis

Since the downloaded RNA-seq expression profile data were converted logarithmically, we first transformed it into a read count generated by normal HTSeq software.13 In addition, the edgeR package14 was used to screen the DEGs between GC tissues samples and normal controls in the expressing data. Next, information on the obtained genes was supplemented by combining the gene annotation information in GENCODE database.15

Identification of co-expression gene of TL1A and GSEA

Pearson’s correlation matrices were performed for the analysis of the correlation between TL1A and other genes. Genes with correlation coefficients >0.5 (P<0.01) were considered to have co-expressed relationship with TL1A. Moreover, according to gene-type information, lncRNA was selected as the key gene for downstream analysis.

To obtain the KEGG pathway associated with TL1A, we conducted GSEA, as previously described,16 based on the background of all the genes. Results with ES >0.5, NES >1.7 and P-value <0.05 were regarded as significant. Gene sets derived from the enrichment analysis were integrated with mRNA from the co-expression analysis to identify overlapped genes as key genes.

Construction of the tl1a-related miRNA-lncRNA/mRNA network

In order to construct the miRNA-lncRNA/mRNA network based on the relationship between miRNA and lncRNA or mRNA, we downloaded the lncRNA transcript sequence, mRNA sequence of 3ʹUTR and mature miRNA sequences from the Ensembl database to predict the relationship between miRNA and lncRNA or mRNA using miRanda.17 Furthermore, the interactions were in accordance with the following criteria: 1) miRNA had a regulatory relationship with TL1A; 2) miRNA simultaneously regulated both mRNA and lncRNA; and 3) miRNA regulated more than 10 mRNAs or lncRNAs. Next, Cytoscape was used to construct and visualize the miRNA-lncRNA/mRNA regulatory network.18

Validation of DEGs based on clinical samples and cell lines

For further validation of our findings, TL1A expression levels were also determined in clinical 10 GC tissues samples and corresponding adjacent normal samples from the Affiliated Hospital of Chengde Medical College, as well as normal human stomach cells GES-1 and a panel of GC cells (AGS, BGC823, MGC803, SGC7901, and HGC27), which were purchased from the Chinese Academy of Sciences, Shanghai, China. This study was approved by the Ethical Committee of Affiliated Hospital of Chengde Medical College, and all patients signed informed consents.

Total RNA from tissues and cells were extracted using Trizol reagent (Invitrogen). First-strand cDNA was generated using Reverse Transcription Reagents (Takara RR047) according to the manufacturer’s protocol. Real-time qPCR was performed in the Real-Time PCR Detection System (Bio-Rad) using SYBR Green (Takara RR820). Primer sequences were as follows: TL1A, F: 5′-TTAGAGCAGACGGAGATAAGCCAAG-3′; R: 5′-CCTAGTTCATGTTCCCAGTGCAGA-3′ and β-actin, F: 5′-AGCGAGCATCCCCCAAAGTT-3′; R: 5′-GGGCACGAAGGCTCATCATT-3′. Relative expression of TL1A normalized to β-actin was calculated with the 2–ΔΔCt method.

Survival analysis

To evaluate the prognostic value of the screened genes, survival analysis was performed, and to improve the reliability of the results with larger samples, we performed a KM plot analysis (http://kmplot.com) using a JetSet best probe set (Only including HGU133 PLUS 2.0 microarray data), based on the datasets, GSE14210, GSE15459, GSE22377, GSE29272, GSE51105, and GSE62254. A P<0.05 was considered statistically significant.

Results

Expression pattern of TL1A in GC

As shown in Figure 1, multidimensional scaling (MDS) analysis of the samples revealed that the tumor samples were concentrated in the middle of the image in compared to the normal controls. Moreover, TL1A expression in GC samples was higher than that in normal controls (Log FC =1.07), and the difference was significant (P=8.90E-07; Table 1).

Figure 1

Multidimensional scaling (MDS) analysis of samples. Compared to normal tissue samples (green circle), tumor samples (red square) were concentrated in the center of the image.

Table 1

The information of TL1A

Categories	Terms
Gene_id	ENSG00000181634.7
Gene_name	TNFSF15
Gene_type	protein_coding
Gene_status	KNOWN
Level	2
Chrome	chr9
Start	1.15E+08
End	1.15E+08
Strand	–
Log FC	1.073987
Log CPM	4.415753
LR	24.15204
P-value	8.90E-07
FDR	8.01E-06

Abbreviations: LR, likelihood ratio tests; FDR, False discovery rate.

Co-expression analysis

A total of 215 genes co-expressed with TL1A were obtained after Pearson’s correlation analysis. As illustrated in Table 2, the top ten genes more closely correlated with TL1A, including caspase 10 (CASP10), iodotyrosine deiodinase (IYD), major vault protein (MVP), g protein-coupled receptor class c group 5 member A (GPRC5A), solute carrier family 45 member 4 (SLC45A4), interferon gamma receptor 1 (IFNGR1), rho guanine nucleotide exchange factor 5 (ARHGEF5), phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit beta (PIK3CB), lymphocyte antigen 75 (LY75), and caspase 8 (CASP8) were identified (all, P<0.05).

Table 2

Analysis of co-expression genes with TL1A (top 10)

Gene ID	Gene name	Gene_type	Cor	P-value	adj. p
ENSG00000003400.13	CASP10	protein_coding	0.586083	6.59E-39	5.13E-35
ENSG00000009765.13	IYD	protein_coding	0.522448	6.98E-30	3.80E-27
ENSG00000013364.17	MVP	protein_coding	0.504588	1.12E-27	3.57E-25
ENSG00000013588.5	GPRC5A	protein_coding	0.562424	2.53E-35	6.65E-32
ENSG00000022567.8	SLC45A4	protein_coding	0.517492	2.94E-29	1.37E-26
ENSG00000027697.11	IFNGR1	protein_coding	0.501095	2.93E-27	8.54E-25
ENSG00000050327.13	ARHGEF5	protein_coding	0.538859	4.99E-32	5.45E-29
ENSG00000051382.7	PIK3CB	protein_coding	0.507103	5.58E-28	1.93E-25
ENSG00000054219.10	LY75	protein_coding	0.560553	4.73E-35	1.14E-31
ENSG00000064012.20	CASP8	protein_coding	0.550115	1.43E-33	2.35E-30

Abbreviations: Cor, correlation coefficient; adj.p, Corrected p-value.

Functional prediction based on GSEA

We performed GSEA to identify the potential function of the screened DEGs. A total of 21 KEGG pathways were enriched (ES >0.5, NES >1.7, P<0.05; Table 3). Furthermore, based on the enriched gene sets and the genes from the co-expression analysis, we identified the intersection genes as the key genes that are relevant to the disease progression of GC. Moreover, the relationships between the intersection genes and the enriched KEGG pathway were visualized as interaction networks (Figure 2).

Table 3

KEGG pathways significantly enriched from GSEA

KEGG	Size	ES	NES	NOM p-val
RIG I like receptor signaling pathway	70	0.629748	1.981377	0
ABC transporters	44	0.614937	1.864473	0.0018622
Cytosolic DNA sensing pathway	55	0.582942	1.83347	0.00566038
Glycosylphosphatidylinositol GPI anchor biosynthesis	25	0.755047	1.795487	0.00782779
Glycerophospholipid metabolism	76	0.585479	1.777131	0.01529637
Regulation of autophagy	35	0.531891	1.707149	0.01651376
JAK-STAT signaling pathway	151	0.523395	1.736979	0.01792829
NOD like receptor signaling pathway	61	0.628086	1.759142	0.01901141
Fc gamma r mediated phagocytosis	96	0.609463	1.74532	0.02294455
N-Glycan biosynthesis	46	0.707322	1.729733	0.02385686
O-Glycan biosynthesis	27	0.626263	1.737538	0.02420857
Sphingolipid metabolism	39	0.624162	1.74533	0.02830189
Leukocyte transendothelial migration	116	0.549586	1.734766	0.02868069
Phosphatidylinositol signaling system	76	0.607689	1.7178	0.02929688
Tight junction	130	0.523988	1.712578	0.03201507
Lysosome	119	0.648934	1.752603	0.03333334
Epithelial cell signaling in Helicobacter pylori infection	67	0.599525	1.721398	0.03578154
Apoptosis	86	0.611604	1.760369	0.03591682
Axon guidance	128	0.545902	1.706259	0.04339623

Figure 2

Relationship between KEGG pathway and enrichment genes from GSEA. Red diamond represents enriched KEGG pathway, while the blue circular node represents the enriched gene of the related pathway.

Construction of miRNA-lncRNA/mRNA network

Based on the above data, the miRNA-lncRNA/mRNA network was constructed using Cytoscape. As illustrated in Figure 3, the miRNA-lncRNA/mRNA network was composed of 7 miRNA nodes, 27 lncRNA nodes, and 21 mRNA nodes.

Figure 3

The miRNA-lncRNA/mRNA network in GC. Blue triangles represent miRNA. Orange diamonds represent lncRNA. Circles represent mRNA.

Identification of survival-related lncRNAs in GC

In order to investigate the genes associated with the survival of GC patients, 876 patients were divided into the high and low groups based on the gene expression median. The results showed that the patients with high TL1A levels had shorter overall survival (OS) than those with low expression (P=2.6e−07; Figure 4A). Similarly, SRC Proto-Oncogenes (SRC) were positively correlated with the OS of GC patients, respectively (P=0.0081; Figure 4B).

Figure 4

Survival analysis using the Kaplan-Meier curve. Patients with high expression of TL1A (A) and SRC (B) had worse OS than those with low expression.

Validation of TL1A in clinical tissue samples and cell lines

To confirm the reliability of the TL1A identified in silico, the TL1A expression profiles were verified in tissue samples and cell lines. TL1A levels in ten GC tissues were all higher than those in adjacent normal samples (Figure 5), which is in accordance with the data above. Similarly, TL1A expression in GC cell lines were all significantly upregulated compared to that in normal cell lines (all, P<0.05).

Figure 5

TL1A expression is determined by qRT-PCR. The TL1A mRNA levels in GC tissues (A) and cell lines (B) were higher than those in normal tissue samples and cell lines. *P<0.05 and **P<0.01.

Discussion

Despite developments in GC treatment, its prognosis remains frustratingly poor. Hence, the identification of risk factors and accurate prognostic prediction in GC patients is still urgent. In this study, we determined the clinical significance of TL1A in GC progression based on TCGA database. The results showed that TL1A expression in GC tissues samples was significantly increased compared to that in normal controls. Moreover, a set of co-expressed genes, enriched genes, and KEGG pathways were obtained. Then, the miRNA-lncRNA/mRNA network was constructed that provided an integrated view of RNA regulation. After a series of bioinformatics analyses, TL1A was identified to be associated with the progression and prognosis of GC.

In fact, as a counter-balancing factor in the maintenance of vascular integrity, TL1A has been found to be abnormally expressed in multiple tumor types.19,20 Moreover, Qin et al21 reveals that TL1A could promote lymphatic metastasis of lung cancer cells via VEGFC up-regulation. These studies suggest that TL1A may play an important role in tumorigenesis, and in accordance with them, our results, based on the data from TCGA database, suggest that TL1A expression levels in GC tissue samples increased significantly when compared to normal controls. Furthermore, the validation analysis showed that TL1A expression was upregulated in other clinical GC tissue samples and cell lines. In addition, co-expression analysis detected 215 genes that were co-expressed with TL1A. Among these co-expressed genes, CASP10 is reported to be associated with the risk of GC,22 and GPRC5A overexpression could predict the advanced poor prognosis in patients with GC.23 All these results suggested that TL1A, along with a series of co-expressed genes are involved in the development of GC; however, the mechanism of function still requires further investigation.

Recent studies suggest TL1A can suppress proliferation and induce apoptosis of endothelial cells through different signaling pathways.24,25 Based on KEGG pathway analysis with GSEA, some cancer-related pathways, such as the phosphatidylinositol, NOD-like receptor (NLR), and JAK-STAT signaling pathways, were determined, indicating that TL1A may play an important role in GC development through different pathways. Chiappini et al26 reported the prognostic implications of the phosphatidylinositol signaling pathway activation in GC. Castaño-Rodríguez et al27 highlighted the relationship between the NLR signaling pathway and gastric carcinogenesis. Additionally, the involvement of the JAK-STAT signaling pathway in gastric carcinogenesis has been reported.28,29 Hence, via pathway analysis, we further verified that TL1A might be an important factor in GC development.

Until now, few studies have investigated the interaction between miRNA and lncRNAs, or miRNAs and mRNAs in GC. In this study, we constructed a miRNA-lncRNA/mRNA network with 7 miRNAs, 27 lncRNAs, and 21 mRNAs, which provided an integrated view of the regulatory crosstalk among these genes associated with TL1A in GC. More importantly, TL1A, along with SRC, were all found to be significantly associated with the survival of GC patients. Specifically, patients with high expression levels of TL1A, or SRC had shorter survival times than those with low expression. In fact, accumulating evidence indicates that TL1A was associated with GC prognosis.11,30 SRC, a nonreceptor tyrosine kinase that acts downstream of receptor tyrosine kinases, is involved in cell proliferation migration, adhesion, and invasion.31 Elsberger et al32 reported that high expression levels of cytoplasmic SRC was associated with a worse prognosis of luminal-type breast cancer, which is in accordance with our studies.

Conclusion

In conclusion, we identified that TL1A was aberrantly expressed in GC; meanwhile, a set of genes co-expressed with TL1A, and pathways were detected. Furthermore, we constructed the miRNA-lncRNA/mRNA network, and identified that TL1A was correlated with OS in GC patients. To our knowledge, this is the first study to evaluate the function of TL1A in the tumorigenesis and progression of GC. However, the clinical value of TL1A for targeted GC therapy and prognosis still needs to be verified.