Innovative signature establishment using lymphangiogenesis-related lncRNA pairs to predict prognosis of hepatocellular carcinoma.

Aims: Hepatocellular carcinoma (HCC) remains a major tumoral burden globally, and its heterogeneity encumbers prognostic prediction. The lymphangiogenesis-related long non-coding RNAs (lrlncRNAs) reported to be implicated in immune response regulation show potential importance in predicting the prognostic and therapeutic outcome. Hence, this study aims to establish a lrlncRNA pairs-based signature not requiring specific expression levels of transcripts, which displays promising clinical practicality and satisfactory predictive capability. Main methods: Transcriptomic and clinical information of the Liver Hepatocellular Carcinoma (LIHC) project retrieved from the TCGA portal were used to find differently expressed lrlncRNA (DElrlncRNA) via analysis performed between lymphangiogenesis-related genes (lr-genes) and lncRNAs(lrlncRNA), and to ultimately construct the signature based on lrlncRNA pairs screened out via Lasso and Cox regression analyses. Akaike information criterion (AIC) values were computed to find the cut-off point optimum for high-risk and low-risk group allocation. The signature then underwent trials in terms of its predictive value for survival, clinicopathological features, immune cells infiltration in tumoral microenvironment, selected checkpoint biomarkers and chemosensitivity. Key findings: A novel lymphangiogenesis-related lncRNA pair signature was established using nine lrlncRNA pairs identified and significantly related to overall survival, clinicopathological features, immune cells infiltration and susceptibility to chemotherapy. Moreover, the signature efficacy was verified in acknowledged clinicopathological subgroups and partially validated by qRT-PCR assay in various human HCC cell lines. Significance: The novel lrlncRNA-pairs based signature was shown to effectively and independently estimate HCC prognosis and help screen patients suitable for anti-tumor immunotherapy and chemotherapy.

Introduction

Despite the rapid advancement of diagnostic and therapeutic techniques, primary liver cancer, of which hepatocellular carcinoma (HCC) constitutes approximately 75%, remains a major constituent to the universal tumor burden (McGlynn et al., 2021). Prevalence of the most meaningful risk factors for HCC at present, hepatitis B and C infection, should drop in years to come owing to vaccination of the newborns (Akinyemiju et al., 2017). However, metabolic risk factors, inclusive of adiposity (Lauby-Secretan et al., 2016), diabetes (Ohkuma et al., 2018), local lymphatic metastasis (Qin and Tang, 2002) and alcoholic abuse (Petrick et al., 2018), are displaying increasing importance and tend to jointly become the major cause of HCC worldwide. This encumbers HCC prognosis prediction due to high heterogeneity of HCC and different risk factors impacting the disease advancement (Torrecilla et al., 2017), thus necessitating the identification of new biomarkers useful for better prognostic prediction and treatment.

Lymphangiogenesis, the process of lymphatic vessel formation, is deeply involved in homeostasis, metabolism and immunity (Suzuki-Inoue et al., 2020). More specifically, tumor-associated lymphangiogenesis plays an essential role in tumor pathogenesis and metastasis through mechanisms like providing niches for tumor stem cells and inhibiting antitumor immune responses (Hu and Luo, 2018). A previous study has revealed that lymphangiogenesis exerted a consequential influence on the survival of HCC patients (Thelen et al., 2009). Some lymphangiogenic genes have been used as biomarkers to predict the prognosis of patients with colorectal liver metastasis after partial hepatectomy (Vellinga et al., 2017). However, prediction based on messenger RNA could suffer from unsatisfactory accuracy due to its inadequate tissue specificity (Deveson et al., 2017). Therefore, it is of necessity to develop new lymphangiogenesis-related biomarkers for better prognostic prediction of HCC.

Long non-coding RNAs (lncRNAs), transcripts whose length is greater than 200 nucleotides, function by regulating gene expression at the post-transcriptional level instead of coding functional proteins (Statello et al., 2021). Accounting for over two-thirds of human transcriptome, lncRNAs play essential roles in various physiological or pathological processes, including the regulation of lymphatic vasculature (Iyer et al., 2015; Md Yusof et al., 2020). LncRNAs are also active participants in tumorigenesis, as evidenced by their roles in liver cancer axis (Zhao and Lawless, 2013). Previous studies have revealed that lncRNAs could interact with genes encoding products to modify the immune microenvironment, thus regulating tumor immune-cell infiltration and aiding the malignant transformation of tumors (Chen et al., 2017).

Previous evidence has demonstrated that lncRNA-based signatures assumed the predictive and prognostic importance in tumors. Three serum lncRNAs have been used to predict the HCC lymph node metastasis status (Ma et al., 2019). Qianhui Xu et al. constructed a 7-immune-associated lncRNA model able to predict prognosis and immunocheckpoint blockade of HCC(Xu et al., 2021). Further, signatures using two-biomarker combinations were proved to be more sensitive and applicable than the ones mentioned, partly because the latter required the specific expression of chosen lncRNAs to be normalized, which is necessary for reducing batch effects among platforms in order to qualify for clinical application (Lv et al., 2020). For example, Ranran Zhou et al. used a ferroptosis-related lncRNA signature consisting of 22 lncRNA pairs to estimate bladder cancer prognostic and immune features (Zhou et al., 2021). However, the number of studies using lncRNA pairs signature for tumoral prognosis prediction is relatively limited. In the present study, we established an innovative prognosis signature using nine lymphangiogenesis-related lncRNA pairs, which efficiently predicted patient survival, tumoral infiltration of immune cells, immunocheckpoint genes expression and chemosensitivity.

Materials and methods

Retrieval of transcriptomic and clinical data

Transcriptomic profiling and clinical information of patients with hepatocellular carcinoma were retrieved from the LIHC project in the TCGA portal2. Samples with a follow-up time less than 30 days or absent survival status were removed. GTF files used for annotation and thus distinguishment between LncRNAs and mRNAs were attained from the ENSEMBL database3. 75 lymphangiogenesis-related genes with relevance scores greater than 1.5 were selected from the GeneCards database4 which was shown in Table 1. Selection of lymphangiogenesis-related genes was justified using Gene Ontology (GO) analysis.

Table 1

Lymphangiogenesis-related genes list.

Gene Symbol	Description	Category	Gifts	GC Id	Relevance score
VEGFC	Vascular Endothelial Growth Factor C	Protein Coding	46	GC04M176683	16.69894028
FLT4	Fms Related Receptor Tyrosine Kinase 4	Protein Coding	50	GC05M180607	16.69141579
VEGFD	Vascular Endothelial Growth Factor D	Protein Coding	33	GC0XM015345	13.53529739
CALCRL	Calcitonin Receptor Like Receptor	Protein Coding	44	GC02M187341	8.252692223
PDPN	Podoplanin	Protein Coding	39	GC01P013583	8.003945351
VEGFA	Vascular Endothelial Growth Factor A	Protein Coding	47	GC06P043770	7.842282772
LYVE1	Lymphatic Vessel Endothelial Hyaluronan Receptor 1	Protein Coding	41	GC11M010753	7.152224541
KDR	Kinase Insert Domain Receptor	Protein Coding	52	GC04M055078	7.094721794
PROX1	Prospero Homeobox 1	Protein Coding	42	GC01P213983	6.581964493
PTPN14	Protein Tyrosine Phosphatase Non-Receptor Type 14	Protein Coding	43	GC01M214348	6.412934303
SOX18	SRY-Box Transcription Factor 18	Protein Coding	39	GC20M064047	4.160199642
PTGS2	Prostaglandin-Endoperoxide Synthase 2	Protein Coding	48	GC01M186640	4.110999107
FLT1	Fms Related Receptor Tyrosine Kinase 1	Protein Coding	49	GC13M028300	3.471794605
CCBE1	Collagen And Calcium Binding EGF Domains 1	Protein Coding	40	GC18M059430	3.466086864
ANGPT2	Angiopoietin 2	Protein Coding	44	GC08M006499	3.383462429
FOXC2	Forkhead Box C2	Protein Coding	44	GC16P086574	3.310676575
NRP2	Neuropilin 2	Protein Coding	43	GC02P205681	3.150856733
FGF2	Fibroblast Growth Factor 2	Protein Coding	46	GC04P122826	2.654223204
PGF	Placental Growth Factor	Protein Coding	42	GC14M074941	2.635710478
TGFB1	Transforming Growth Factor Beta 1	Protein Coding	50	GC19M041301	2.578747272
CCR7	C–C Motif Chemokine Receptor 7	Protein Coding	44	GC17M040556	2.554449558
HIF1A	Hypoxia Inducible Factor 1 Subunit Alpha	Protein Coding	46	GC14P061695	2.478032827
HMGB1	High Mobility Group Box 1	Protein Coding	44	GC13M030456	2.41782546
CXCR4	C-X-C Motif Chemokine Receptor 4	Protein Coding	51	GC02M136114	2.304786205
POSTN	Periostin	Protein Coding	42	GC13M037562	2.292275667
CCL21	C–C Motif Chemokine Ligand 21	Protein Coding	41	GC09M034709	2.277864933
VASH1	Vasohibin 1	Protein Coding	36	GC14P076761	2.217816591
SHH	Sonic Hedgehog Signaling Molecule	Protein Coding	48	GC07M155799	2.203608274
STAT3	Signal Transducer And Activator Of Transcription 3	Protein Coding	51	GC17M042313	2.19051075
ANGPT1	Angiopoietin 1	Protein Coding	45	GC08M107246	2.184707403
NRP1	Neuropilin 1	Protein Coding	45	GC10M033177	2.166389465
ERBB2	Erb-B2 Receptor Tyrosine Kinase 2	Protein Coding	52	GC17P039687	2.150751829
NR2F2	Nuclear Receptor Subfamily 2 Group F Member 2	Protein Coding	48	GC15P096325	2.107103348
CXCL12	C-X-C Motif Chemokine Ligand 12	Protein Coding	44	GC10M044294	2.088619471
VEGFB	Vascular Endothelial Growth Factor B	Protein Coding	43	GC11P064234	2.052349567
HPSE	Heparanase	Protein Coding	44	GC04M083292	2.029333591
TEK	TEK Receptor Tyrosine Kinase	Protein Coding	49	GC09P027109	2.020673752
HGF	Hepatocyte Growth Factor	Protein Coding	50	GC07M081699	2.000943422
MET	MET Proto-Oncogene, Receptor Tyrosine Kinase	Protein Coding	52	GC07P116672	1.986150503
CEACAM1	CEA Cell Adhesion Molecule 1	Protein Coding	41	GC19M042507	1.980372667
NFKB1	Nuclear Factor Kappa B Subunit 1	Protein Coding	51	GC04P102501	1.93708396
NOS2	Nitric Oxide Synthase 2	Protein Coding	49	GC17M027756	1.935756326
S1PR1	Sphingosine-1-Phosphate Receptor 1	Protein Coding	44	GC01P101236	1.911855936
FOXC1	Forkhead Box C1	Protein Coding	42	GC06P001610	1.900332332
PDGFB	Platelet Derived Growth Factor Subunit B	Protein Coding	48	GC22M051474	1.89337194
CLEC14A	C-Type Lectin Domain Containing 14A	Protein Coding	33	GC14M038254	1.842555404
SMAD4	SMAD Family Member 4	Protein Coding	49	GC18P051028	1.834085941
IL17A	Interleukin 17A	Protein Coding	41	GC06P052186	1.829650402
ITGA4	Integrin Subunit Alpha 4	Protein Coding	48	GC02P181456	1.812556744
IL7R	Interleukin 7 Receptor	Protein Coding	45	GC05P035852	1.789921403
TNF	Tumor Necrosis Factor	Protein Coding	49	GC06P061170	1.782137632
SIX1	SIX Homeobox 1	Protein Coding	43	GC14M060643	1.782137632
MMP9	Matrix Metallopeptidase 9	Protein Coding	52	GC20P046008	1.76113379
SMARCA4	SWI/SNF Related, Matrix Associated, Actin Dependent Regulator Of Chromatin, Subfamily A, Member 4	Protein Coding	48	GC19P010932	1.741836548
IL6	Interleukin 6	Protein Coding	48	GC07P022725	1.740508795
MIR27B	MicroRNA 27b	RNA Gene	21	GC09P095097	1.720307231
MIR9-1	MicroRNA 9-1	RNA Gene	21	GC01M156420	1.703069806
PECAM1	Platelet And Endothelial Cell Adhesion Molecule 1	Protein Coding	40	GC17M064319	1.695909262
MAPK14	Mitogen-Activated Protein Kinase 14	Protein Coding	50	GC06P061307	1.67373991
MCAM	Melanoma Cell Adhesion Molecule	Protein Coding	39	GC11M119308	1.67373991
TYMP	Thymidine Phosphorylase	Protein Coding	45	GC22M050525	1.671724916
EDN1	Endothelin 1	Protein Coding	46	GC06P012256	1.657824516
ITGB1	Integrin Subunit Beta 1	Protein Coding	49	GC10M032899	1.657404423
PDGFA	Platelet Derived Growth Factor Subunit A	Protein Coding	42	GC07M000497	1.648259282
IL24	Interleukin 24	Protein Coding	41	GC01P206897	1.591114283
TIAM1	TIAM Rac1 Associated GEF 1	Protein Coding	44	GC21M031118	1.558186173
ECM1	Extracellular Matrix Protein 1	Protein Coding	43	GC01P150508	1.558186173
LIMS1	LIM Zinc Finger Domain Containing 1	Protein Coding	40	GC02P108534	1.558186173
NES	Nestin	Protein Coding	39	GC01M156668	1.558186173
CDKN2B-AS1	CDKN2B Antisense RNA 1	RNA Gene	21	GC09P021994	1.558186173
ADM	Adrenomedullin	Protein Coding	44	GC11P010304	1.55810535
ITGA9	Integrin Subunit Alpha 9	Protein Coding	42	GC03P037468	1.546670914
MMP2	Matrix Metallopeptidase 2	Protein Coding	52	GC16P055390	1.536451578
NFATC1	Nuclear Factor Of Activated T Cells 1	Protein Coding	47	GC18P079395	1.511538029
TIE1	Tyrosine Kinase With Immunoglobulin Like And EGF Like Domains 1	Protein Coding	42	GC01P043300	1.511538029

Identification of lymphangiogenesis-related long noncoding ribonucleic acids pairs

Correlation analysis conducted between lymphangiogenesis-related genes and the entirety of lncRNAs was employed to screen lrlncRNAs. The standard for screening was defined as correlation coefficient greater than 0.4 and p-value less than 0.001. Subsequently, differential expression analysis was conducted among lrlncRNAs to identify differently expressed lymphangiogenesis-related lncRNAs (DElrlncRNAs) with the criterion of FDR less than 0.05 and log2 |fold change (FC)| greater than 2. To identify DElrlncRNA pairs, DElrlncRNAs were cyclically paired in a lncRNAa-lncRNAb pattern to form a zero-or-one matrix, in which if the expression of lncRNAa is higher than that of lncRNAb, value 1, otherwise value 0, would be yielded. The established zero-or-one matrix was subjected to further screening until DElrlncRNA pairs whose proportion of being value 1 or 0 was less than 20% or more than 80% were removed. The DElrlncRNA pairs removed are regarded as unnecessary for further analysis in that their expression ratios are considered the same in all samples (Hong et al., 2020).

Construction and verification of the prognostic signature using differently expressed lymphangiogenesis-related long noncoding RNA pairs

Uni-Cox regression was performed for DElrlncRNA pairs left in the matrix with a p-value less than 0.01 to pick out the ones with prognostic significance. LASSO regression analysis was adopted to avert overfitting. Multivariate Cox regression was conducted later to establish the prognosis signature, which was used to calculate risk scores for clinical samples with the formula: risk score = , where n means the quantity of DElrlncRNA pairs within the prognostic signature and val(i) and coef(i) represent the value yielded in the matrix and the regression coefficient, respectively. To evaluate the prognosis signature, ROC curves were plotted with corresponding AUCs calculated and comparisons made between this signature and other clinical variables. Patients were allocated into high- or low-risk group as per the cut-off point generated from the Akaike information criterion (AIC) values of 1-year ROC curve. Survival differences between the two risk groups were compared using Kaplan-Meier method and logrank test. Subsequently, univariate and multivariate analyses were performed to discern independently prognostic predictors of HCC patients, and more importantly, whether the risk score serves as one. The relation between the prognosis signature and other clinicopathological features was investigated using chi-square tests. Differences in risk scores among subgroups with these clinical characteristics were shown as box plots via Wilcoxon signed-rank test. A nomogram model was constructed using the two independently prognostic predictors unanimously identified by the univariate and multivariate analyses and validated by calibration graphs of 1-/2-/3-year comparing the actual survival probability of HCC patients and the one yielded by the nomogram (Wan et al., 2017).

Biofunction and pathways exploration

A coexpression network between lymphangiogenesis-related lncRNAs in the signature and their interactive mRNA previously identified was built using Cytoscape (version 3.9.0). The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted with R packages org. Hs.eg.db (version 3.13.0) and GOplot and visualized by ggplot2, under the criterion of p-/q-values<0.05.

Tumoral infiltration of immune cells analysis

The correlation between the risk score and the tumoral infiltration of immune cells estimated by well-received methods at present inclusive of XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT−ABS and CIBERSORT was investigated for later visualization in the form of a lollipop diagram via Spearman correlation test, whose significance threshold was specified as p-value less than 0.05. Wilcoxon signed-rank test was employed to analyze, and boxplots were plotted to display the infiltration differences of immune cells between the two risk groups.

Immune checkpoint genes expression

For a better understanding of the expression differences in immunocheckpoint genes between the high-risk and low-risk groups, several immunocheckpoint genes were chosen for analysis, including CD47, CD276, LAG3, CTLA4, PDCD1 and HAVCR2; the results are visualized as violin plots by the ggpubr R package.

Chemosensitivity evaluation

To examine the prognosis signature under clinical applications, half-maximal inhibitory concentration (IC50), which indicates the concentration of drug needed to inhibit tumor cells by 50%, was calculated for five chemotherapeutic drugs reportedly useful for patients with hepatocellular carcinoma, including doxorubicin, gemcitabine, mitomycin C and sorafenib. Boxplots were drawn to display the contrasts in the IC50s between high-risk and low-risk groups computed by Wilcoxon signed-rank test via pRRophetic R package.

Cell lines and cell culture

Human normal liver cells (LO2) and human HCC cell lines (HepG2, SK-HEP-1 and Huh7) were purchased from the Cell Bank (Cell Institute, Sinica Academia Shanghai, Shanghai, China) and validated by short tandem repeat (STR) profiling. Cells were either cultured in RPMI 1640 (Gibco, USA) or Dulbecco's Modified Eagle's Media (DMEM, Biological Industries, Israel) supplemented with 10% fetal bovine serum (FBS, Gibco, USA) and 1% penicillin-streptomycin and incubated at 37 °C in 5% CO2.

Verification of LncRNA pairs by quantitative real-time polymerase chain reaction

The EZ-press RNA Purification Kit (EZBioscience, Shanghai, China) was used to extract total RNA from the human normal and tumor cells following the manufacturer's instructions. The purity and concentration of RNA extracted were measured using IMPLEN N60 Touch (IMPLEN, Germany). HiScript II Q RT SuperMix for qPCR (+gDNA wiper) (Vazyme, Nanjing, China) was used for reverse transcription and ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) for qRT-PCR reaction under the LightCycler96 System (Roche, Germany). LncRNA pair expression ratios were calculated using the relative 2 –ΔΔCt method with internal control of β-actin. The differences in the ratios in liver normal and tumor cells were compared by t-tests and visualized in bar graphs drawn with GraphPad Prism (version 9.0). See Table 2 for the list of qRT-PCR primer sequences.

Table 2

QRT-PCR primer sequences.

Species	Gene name	Primer sequence (5→3′)
Homo sapiens	AC068506.1	Forward	TCCCATCTCCCACTATTC
		Reverse	AAGGCACATACAAGAAAGC
Homo sapiens	LENG8−AS1	Forward	AGCACGGACTCTGATACAA
		Reverse	TCAGCCAGTTCTCCCTAAT
Homo sapiens	AC006042.1	Forward	TACTTTTACCCTTGAGCA
		Reverse	GAACATCTACAATGAGCC
Homo sapiens	AL355488.1	Forward	AGCACCTTGGTTCTGATGT
		Reverse	CCTGGCTATGGCACTTACT
Homo sapiens	ACTB	Forward	CATGTACGTTGCTATCCAGGC
		Reverse	CTCCTTAATGTCACGCACGAT

Results

Differently expressed lymphangiogenesis -related long noncoding RNAs

The flow chart of this study was shown as Figure 1. Transcriptomic profiles and clinical information of 377 LIHC patients were retrieved from TCGA portal. After removing the patients whose follow-up time was less than 30 days, 349 patients were included for further analysis. Lymphangiogenesis-related genes were attained from the GeneCards database and subjected to Gene Ontology analysis, in which GO terms justified the gene selection in this fashion, before and after being used for coexpression analysis to screen out lrlncRNAs (Figure 2C). Under the criterion of FDR <0.05 and log2 |fold change (FC)| >2, a total of 32 lrlncRNAs were found, all of which were upregulated (Figure 2A, B).

Figure 1

Flow chart of the study.

Figure 2

DElrlncRNAs identification; Identifying differently expressed lymphangiogenesis-related lncRNAs (DElrlncRNAs) using the LIHC dataset from TCGA portal, as shown in the heatmap (A) and the volcano plot (B) (C) GO terms indicated the selected genes were lymphangiogenesis-related.

Establishment and validation of lymphangiogenesis-related long noncoding RNA pairs and prognosis signature

Using the zero-or-one matrix that singly and cyclically paired DElrlncRNAs, 346 valid DElrlncRNA pairs were constructed, among which 16 pairs were extracted after LASSO regression analysis following a univariate Cox regression analysis (Figure 3A, B). Subsequently, a Cox proportional hazard model was constructed stepwise with 9 of the 16 pairs; results were shown as Figure 3C, D. To elucidate the biological functions and pathways relevant to the established signature, lncRNAs within the signature and their coexpressed mRNAs were used to build the network (Figure 4A). The GO and KEGG analyses therefore performed indicated close functional connection of the signature with chemotaxis and endothelial cell proliferation (Figure 4B), and participation in Rap1, Ras, MAPK, PI3K−Akt, Calcium and AGE−RAGE signaling pathway (Figure 4C). Further, AUC values were calculated for the ROC curves of all the 9 DElrlncRNA pairs, the maximum of which was 0.825 and was used for calculation of the AIC value to yield the optimal cut-off point (Figure 4D). ROC curves of 1-/2-/3-year were plotted, and other clinical features were also compared to assess the optimality of the signature (Figure 4E, F). The AUC of 1-year ROC curve was much greater than that of other clinical characteristics, endorsing the clinical significance the signature possesses based on DElrlncRNA pairs (Figure 4F). Using the cut-off point 1.580, patients were allocated to different risk groups depending upon their risk scores (Figure 5A). It can be seen from the scatterplot and boxplot that patients in the low-risk group had a higher chance of survival (Figure 5B). The conclusion was also supported by the survival curve showing the survival probability of the high-risk patients was significantly lower than that of the low-risk (Figure 5C, p < 0.001). Subgroup survival analyses performed to reduce bias rendered similar results (Figure 6). To explore the relation between the risk score and other clinical features, Chi-square tests were conducted among such subgroups of clinical features as gender, age, clinical stage, T/M/N stage and Grade. The results were displayed in a heatmap presenting the extent to which these clinical features were related to the signature. Of all, clinical stage, Grade and T stage were the ones significantly associated (Figure 7A). Wilcoxon signed-rank tests performed for the same purpose rendered similar results except for M stage (Figure 7B-H), which may be attributed to insufficient cases of M1. To validate the signature without other clinical characteristics affecting the outcome, univariate (Figure 7I) and multivariate (Figure 7J) Cox regression analyses were adopted to discern independently prognostic predictors of LIHC patients and the results from both analyses agreed that risk score (Univariate: p < 0.001, HR = 1.208, CI(1.158–1.259); Multivariate: p < 0.001, HR = 1.190, CI(1.139–1.244)) and stage (Univariate: p < 0.001, HR = 1.808, CI(1.463–2.234); Multivariate: p < 0.001, HR = 1.652, CI(1.324–2.061)) were independent predictors of prognosis. Therefore, the risk score calculated using the signature works decently as recognized clinical predictors. To ameliorate the prediction of LIHC patients' survival, another independently prognostic predictor stage was combined with the risk score to form a nomogram model (Figure 8A). For instance, the estimated 3-year survival rate of a stage III LIHC patient with a risk score of 14 is less than 5%. Calibration graphs of 1-/2-/3-year comparing the actual survival probability of the patients and the one predicted by the nomogram revealed that differences between them were marginal (Figure 8B-D), suggesting the applicability of the nomogram model.

Figure 3

Establishment of prognosis signature (A) 30 lncRNA pairs were analyzed by LASSO regression, in which lncRNA pairs were eventually removed from the model as the penalty (lambda) increased (B) The tuning parameter selection of the LASSO analysis, in which 16 lncRNA pairs were left (−4 < lambda.min<−3.7) (C) The univariate Cox regression analysis of the 9 significant DElrlncRNA pairs used to construct the signature (D) The forest map of the 9 DElrlncRNA pairs screened out by the Cox proportional hazard model, which was also used to select pairs for qRT-PCR validation.

Figure 4

Biofunctions, pathways and predictivity of the signature (A)LncRNA and mRNA coexpression network of the signature (B) Biological Functions identified using GO analysis in the signature (C) Pathways associated with the signature as found out by KEGG analysis (D) The optimal cut-off point was calculated to allocate patients into two different risk groups using the AIC value (E) The AUCs for 1-, 2- and 3-year ROC curves were 0.825, 0.764, and 0.741, respectively (F) The risk score predicted with best efficiency comparing with other clinical features for 1-year survival.

Figure 5

Signature survival prediction Risk scores (A) and survival status (B) of every case were shown and compared (C) The Kaplan-Meier plot showing the significantly slimmer chance of survival for the high-risk group. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001, also applicable for the following figures.

Figure 6

Survival analysis for subgroups of gender (A, B), age (C, D), grade (E, F), stage (G, H), N stage (I, J), M stage (K, L) and T stage (M–O). P-values in all subgroups indicated statistical significance.

Figure 7

Association of the signature with clinicopathological features (A)The heatmap showed that grade, clinical stage and T stage were significantly related to the risk score (B–H) Boxplots using Wilcoxon signed-rank tests agreed that T stage (B), stage (C) and grade (D) significantly correlated with the risk score, while M stage (E), age (F) and gender (G)did not, with the exception of N stage (H). Univariate (I) and multivariate (J) Cox regression analyses discerning independently prognostic predictors.

Figure 8

Using nomogram to predict patients' survival (A)The nomogram model using the risk score and clinical stage to predict 1-/2-/3-year survival rates of LIHC patients (B–D) Calibration graphs indicating the predicted survival rates of 1- (B), 2- (C) and 3-year (D) nomogram model were comparable to the actual ones.

Exploration of relation between risk score and tumoral infiltration of immune cells

Given the involvement of lncRNA in the regulation of tumoral infiltration of immune cells, the relation between the lncRNA-based risk score and immune cells infiltration was explored using Wilcoxon signed-rank tests and Spearman correlation. A detailed list of immune cell types with significant Spearman correlation coefficients was presented in Figure 9A. It was revealed that high risk correlated with higher tumoral infiltration of macrophages (Figure 9B), Th2 cells (Figure 9C), myeloid dendritic cells (Figure 9D), Treg cells (Figure 9E) and neutrophils (Figure 9F), and with lower tumoral infiltration of CD8+ naïve cell (Figure 9G), hematopoietic stem cells (Figure 9H), endothelial cells (Figure 9I) and central memory T cells (Figure 9J).

Figure 9

Association of the signature with tumoral infiltration of immune cells (A) The correlation of the risk score with many types of tumor-infiltrating cells (B) High risk correlated with higher tumoral infiltration of macrophages (B), Th2 cells (C), myeloid dendritic cells (D), Treg cells (E) and neutrophils (F) and lower infiltration of CD8+ naïve cell (G), hematopoietic stem cells (H), endothelial cells (I) and central memory T cells (J).

Risk-related expression of immunocheckpoint genes

To find out if the risk score could be used to predict immune checkpoint blockage therapy, the expression of immunocheckpoint genes in the high-risk and low-risk groups were visualized comparatively in violin plots. It was shown that the expression of HAVCR2 (Figure 10A, p < 0.001), CD47 (Figure 10B, p < 0.001) and CD276 (Figure 10C, p < 0.001) were significantly higher in high-risk group while the expression difference of the rest of the genes analyzed (Figure 10B-F, p > 0.05) showed no statistical significance.

Figure 10

Association of the signature with immunocheckpoint genes; Expression of HAVCR2 (A), CD47(B) and CD276(C) was significantly higher in the high-risk group, while difference in LAG3 (D), PDCD1 (E) and CTLA4 (F) expression displayed no statistical significance between the groups.

Chemotherapeutic prediction using risk score

IC50 values evaluating the chemosensitivity in LIHC patients were computed and the differences between the high-risk and low-risk groups were analyzed using Wilcoxon signed-rank test. As shown in the boxplots, the half inhibitory concentration of gemcitabine (Figure 11A, p < 0.001), doxorubicin (Figure 11B, p < 0.01), mitomycin C (Figure 11C, p < 0.001) and sorafenib (Figure 11D, p < 0.001) was significantly lower in the high-risk group, implying the promise of this signature to estimate the chemotherapeutic outcome.

Figure 11

Association of the signature with chemosensitivity; IC50 of gemcitabine (A), doxorubicin (B), mitomycin C (C) and sorafenib (D) were significantly lower in the high-risk group.

Validating expression-ratios of RNA pairs by quantitative real-time polymerase chain reaction

To validate the ratios of the lncRNA pairs expression, qRT-PCR was performed in normal and tumoral liver cell lines. As shown in Figure 12, the expression ratio of AC068506.1| LENG8-AS1 was significantly elevated, while the ratio of AC006042.1| AL355488.1 was significantly decreased in the high-risk group, according with their hazard ratios in the signature, which indicated that the lncRNA pairs are worthy of further investigation.

Figure 12

Verification of expression ratios in lymphangiogenesis-related lncRNA pairs; qRT-PCR results for the expression ratios of AC006042.1| AL355488.1 (A) and AC068506.1| LENG8-AS1 (B).

Discussion

As a major cause of tumor-related mortality globally, HCC can be highly metastatic and recurrent, which restricts patients' long-term survival (Vogel et al., 2019). Lymphangiogenesis, which describes the growth of new lymphatic vessels, has been shown to relate to metastases and unsatisfactory prognosis in a variety of human tumors, including melanoma, prostate and breast cancers (Rinderknecht and Detmar, 2008). More importantly, recent evidence has revealed that lymphangiogenesis facilitates metastasis in HCC (Yu et al., 2017). As lymphangiogenesis is rarely observed in healthy adults, therapies targeting lymph vessel formation should have the advantage of not intervening normal physiology (Mumprecht and Detmar, 2009). Anti-lymphangiogenic strategies have been developed over the decade to hinder lymphatic metastasis and currently proceed to the stage of clinical trials (Dieterich and Detmar, 2016). There is sufficient experimental evidence that drugs blocking the lymphangiogenic axis reduce tumor metastasis, lymphatically and distantly (Burton et al., 2008; Caunt et al., 2008). Together with the recent discovery that lymphangiogenesis regulates specific immune responses, it is tempting to develop lymphangiogenesis-related biomarkers as latent diagnostic and therapeutic targets for patients with HCC. In recent years, signatures based on specific expression of certain transcripts were proposed to predict the prognosis of malignancy, which required complicated calibration before clinical application (Wu et al., 2021; Xia et al., 2021; Xu et al., 2021). In the present study, we established a lymphangiogenesis prognostic signature for HCC patients by taking advantage of the relative expression of lncRNA pairs, thus allowing better practicability for clinical use.

First, differently expressed lrlncRNAs were sifted using correlation analysis conducted between lr-genes and the whole of lncRNAs with the data retrieved from TCGA portal. Next, valid lrlncRNA pairs were sifted using the zero-or-one matrix and the ones with prognostic significance were screened out and the prognostic signature was established after a series of computations. To test the signature, ROC curves of 1-/2-/3-year were plotted and compared to that of other clinical features such as gender, age and stage. The patients were allocated into high-risk or low-risk group as per the cut-off point computed using the AIC value. Subsequently, the relationship was investigated between the risk score and survival, clinical features, tumoral infiltration of immune cells, expression of immunocheckpoint inhibitor genes and chemosensitivity.

Some of the differently expressed lrlncRNAs included in the signature have already been shown to be important players in HCC. Recent evidence indicated that LINC00205 promoted proliferation of HCC cells by targeting miR-122–5p or miR-26a-5p (Zhang et al., 2019; Cheng et al., 2021). LINC00665 was reported to increase malignancy of HCC through the activation of NF-κB signaling (Ding et al., 2020). LncRNA MYLK-AS1 was also found to facilitate tumor progression of HCC through miR-424–5p/E2F7 axis or EGFR/HER2-ERK1/2 signaling pathway (Liu et al., 2020; Teng et al., 2020). Moreover, expression ratios of AC068506.1| LENG8-AS1 and AC006042.1| AL355488.1 were validated by qRT-PCR in various HCC cell lines, which along with the above studies suggests that the lrlncRNAs within the established signature can be worthy of further investigation.

Tumor microenvironment embraces a broad spectrum of intricate interactions between immune cells, tumor cells and stroma, in which lymphangiogenesis also plays important roles in regulating antitumor immunity (Marin-Acevedo et al., 2018; Garnier et al., 2019). Immune cell activation and infiltration in HCC affect response to anti-tumor blockade and relate to prognosis and therapeutic efficacy (Kurebayashi et al., 2018). To investigate the relation between tumoral infiltration of immune cells and risk score, the following acknowledged methods were used, including TIMER (Li et al., 2017), XCELL (Aran et al., 2017), QUANTISEQ (Finotello et al., 2019), CIBERSORT-A (Tamminga et al., 2020), CIBERSORT (Newman et al., 2015), MCPcounter (Dienstmann et al., 2019) and EPIC (van Veldhoven et al., 2011). In particular, the result revealed that high risk correlated with higher tumoral infiltration of macrophages, Th2 cells, myeloid dendritic cells, Treg cells and neutrophils and lower infiltration with CD8+ naïve cell, hematopoietic stem cells, endothelial cells and central memory T cells. A study recently demonstrated the possibility of predicting the therapeutic benefits of immunotherapy and chemotherapy using the immunogenomic analysis-derived immune scores (Dai et al., 2020). The lrlncRNA signature proposed that high risk was linked with chemosensitivity to such therapeutics as cisplatin, doxorubicin, gemcitabine, mitomycin C and sorafenib, and with high expression of CD47, CD276 and HAVCR2. CD276 molecule, an immune checkpoint also known as B7–H3, inhibits anti-tumor immunity and promotes progression. In the tumor microenvironment, B7–H3 was capable of inhibiting Th1 activation and promoting Th2 differentiation (Feng et al., 2021). The shift from Th1 to Th2 was reported to promote cancer progression (Kumar et al., 2017). The effect of regulatory T cells on anti-tumor immune response suppression has been established (Tanaka and Sakaguchi, 2017). Dendritic cells exposed to regulatory T cells upregulate B7–H3 expression, producing an inhibitory phenotype (Mahnke et al., 2007). A plausible hypothesis was thus proposed that immunotherapy targeting B7–H3 allows the rebalance of Th1/Th2 and reversal of DC inhibitory phenotype. The results in the study were in tone with previous evidence and together these findings suggest that the signature could potentialize clinical immunotherapy and chemotherapy guiding for patients with HCC.

However, we acknowledged as well some weaknesses and limitations in the present article. The raw data for initial analyses were relatively limited since it was merely retrieved from the TCGA portal due to the unavailability of datasets with simultaneous inclusion of lncRNA expression, clinicopathological features, and survival endpoints for HCC patients. Moreover, owing to the expression difference in samples, which might bring unsoundness to the ultimate model, external data validation was required, regardless of the zero-or-one matrix constructed to diminish sample errors sourcing and methods used for optimality testing. Subsequent molecular biological experiments are necessary and under consideration to further investigate the roles of DElrlncRNAs in HCC advancement and the underlying mechanisms. To improve the prognosis-predictive value of the signature, clinical cases would be preferred, the recruitment and analysis of which are time-consuming, even though some of the lncRNA pairs were validated in various human HCC cell lines.

In conclusion, the present study demonstrated that a novel lrlncRNA-pairs based signature without test platforms limitations and not requiring specific expression levels of selected transcripts potentialize prognosis prediction of HCC and may help screen patients suitable for anti-tumor immunotherapy and chemotherapy.

Declarations

Author contribution statement

Jincheng Cao: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Yanni Xu: Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Xiaodi Liu: Conceived and designed the experiments; Analyzed and interpreted the data; Wrote the paper.

Yan Cai: Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data.

Baoming Luo: Conceived and designed the experiments; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

Prof Baoming Luo was supported by National Natural Science Foundation of China [82171944 & 81873899], Natural Science Foundation of Guangdong Province [2021A1515012611].

Data availability statement

Data associated with this study has been deposited at The Cancer Genome Atlas; Ensembl genome database; GeneCards database.

Declaration of interest’s statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.