Deciphering the role of ELAVL1: Insights from pan-cancer multiomics analyses with emphasis on nasopharyngeal carcinoma

Pan-cancer data collection and processing

The UCSC Xena platform was utilized to access the TCGA and GTEx databases, focusing on pan-cancer ELAVL1 expression levels and associated clinical characteristics. Abbreviations of the 33 cancer types encompassed in this study were referred to Supplementary TableS1. We obtained the pan-cancer scRNA-seq datasets from the TISCH database.^[27] Pan-cancer proteomics and ST datasets were summarized and obtained from Sparkle database (https://grswsci.top/) and the corresponding ID for each ST sample were summarized in Supplementary Table S2. NPC datasets were downloaded from the GEO database. ^{[28, 29, 30, 31, 32]} The probes were mapped utilizing the “AnnoProbe” R package, and averaging multiple probes were calculated using the “limma” R package when necessary.^[33]

Differential expression and localization analyses

The HPA database was utilized to assess the ELAVL1 expression level in cell lines and tissues. Data from the TCGA and GTEx databases were obtained to compare differences between tumors and normal tissues at mRNA levels, and pan-cancer proteomics datasets were used to compare differences between tumors and normal tissues at protein levels. The 3D configuration of ELAVL1 was obtained from the PDB database (https://www.rcsb.org/) and the location of ELAVL1 was acquired from the UniProt database (https://www.uniprot.org/).

Diagnostic and prognostic analyses

The ROC curves for the cancers of interest were created using the “pROC” R package. Survival curves were conducted using “survival” and “survminer” R packages. Univariate Cox regression analyses were executed using the “survival” and “forestplot” R packages to assess the prognostic significance of ELAVL1 expression as a predictor of OS, DSS, PFI, and DFI.

Genomic alteration analyses

Pan-cancer analyses of genomic mutation frequencies, amplifications, and deep deletions were performed utilizing the Cancer Type Summary module available in the cBioPortal database.^[34] The processed SNV data were analyzed using the “maftools” R package, in order to elucidate the mutational landscape of ELAVL1 across pan-cancer datasets.^[35] The processed CNV data were obtained from GSCA database and visualized using “ggplot2” R package.^[36]

Modification analyses

PTMs was obtained from the PhosphoSite database.^[37] Differential expression analyses of phosphorylation and methylation modifications were downloaded from the UALCAN database.^[38]

Alternative splicing analyses

The investigation of clinically significant AS of ELAVL1 was conducted utilizing the OncoSplicing server and its ClinicalAS tool.^[39] We identify pertinent ELAVL1 AS events within the SplAdder function.

Functional enrichment and immunological interaction analyses

Based on ELAVL1 mRNA expression levels, the top quarter of samples were defined as the high expression group, and the bottom quarter were considered as the low expression group in each tumor type. We utilized the “limma” R package to calculate the log2 fold change (log2FC) for each gene, and subsequently, all genes were ranked based on their log2FC values.^[33] GSEA was executed using the “clusterProfiler” R package and the hallmark gene set.^[40] Genome-wide screening CRISPR were downloaded from DepMap database,^[41] and candidate genes were calculated using CERES algorithm. A negative score means cell growth inhibition and/or death after gene knockout.

For immunological interaction analyses, we applied “immunedeconv” R package to calculate the infiltration levels from the TCGA pan-cancer datasets,^[42] and immune-related genes were collected from a previous study.^[43] The verification of genome status and immune subtypes were originated from another previous study.^[44] Pan-cancer immunotherapy cohorts were obtained from the BEST database.^[45]

Cell lines and cultures

Human NPC cell lines (CNE1, CNE2, HNE1, C666-1, SUNE-1, and HONE-1) were maintained in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS). N₂-Tert was cultured in keratinocyte serum-free medium (Invitrogen), with supplementation of bovine pituitary extract. The cell lines mentioned above were authenticated and kindly provided by Prof. Wei Xiong (Central South University, China). Radiotherapy-sensitive and resistant HK1 cell lines were kindly provided by Prof. Jian Zhang (Southern Medical University, China).^{[46, 47]}

Western blotting

Cell lysates were prepared using RIPA lysis buffer to extract proteins, which were then separated by SDS-PAGE (Beyotime) and transferred onto PVDF membranes (Millipore). Following 1 hour of blocking in nonfat milk, membranes were incubated with a primary antibody overnight at 4°C, followed by incubation with secondary antibodies at room temperature for 1 hour. The following primary antibodies were used: anit-ELAVL1 (Rabbit, 11910–1-AP, Proteintech), and anti-α-Tubulin (Rabbit, 66031-1-Ig, Proteintech).

Transient transfection and stable cell line construction

For siRNA-mediated RNA interference, siRNAs specifically targeting ELAVL1, as well as non-targeting control, were synthesized and procured from GenePharma (Suzhou, China). The siRNAs were transfected into the cells using Lipofectamine 3000 (Invitrogen), following the protocol provided by the manufacturer. For stable cell line construction, NPC cells were transfected with the shRNA‐ELAVL1 or control viruses. Targeted sequences for siRNA and shRNA are provided in Supplementary TableS3 and Supplementary TableS4.

CCK-8 assays

A total of 2, 000 cells were incubated in 96‐well plates. Subsequently, the CCK-8 solution was introduced into each well and incubated for a duration of 2 h for evaluating cell viability. OD450 values were measured over 5 days.

Transwell assays

A total of 4 × 10⁴ HONE-1 cells and 6 × 10⁴ SUNE-1 cells underwent enzymatic digestion and were subsequently resuspended. Cells were introduced into the upper chambers devoid of FBS, whereas the lower cross-pore compartment contained a medium supplemented with 20% FBS. After a 20-hour incubation, the migrated NPC cells were fixed with methanol and stained with 0.1% crystal violet.

Wound healing assays

The transfected NPC cells were cultured in 6-well plates, and wounds were created using a 100 μL pipette tip. Microscopic images were taken at 0 and 24 hours. ImageJ software was utilized to quantify the scratch area.

Colony formation assays

A total of 1, 000 NPC cells were seeded into each well of 6-well plates using a medium supplemented with 10% FBS. After 2 weeks, the colonies were fixed using methanol for 20 minutes and stained with 0.1% crystal violet for 20 minutes.

Animal experiments

Animal experiments were performed in accordance with the guidelines of the Institutional Animal Care and Use Ethics Committee of SYSUCC (L102012024120B). Female BALB/c nude mice aged 4 weeks were purchased from Zhuhai BesTest Bio‐Tech Co. Ltd. and housed at the Animal Facility of SYSUCC under controlled conditions.

For the popliteal lymph node metastatic model, ELAVL1-KD or control HONE-1 cells (1 × 10⁶ cells) in PBS were subcutaneously injected into the footpads of mice. After one month, the mice were humanely euthanized. The primary footpad along with lymph node metastatic tumors were surgically excised.

For the lung metastasis model, ELAVL1-KD or control HONE-1 cells (1 × 10⁶ cells) were administered via tail vein injection. After one month, the mice were humanely euthanized, and their lungs were surgically excised.

For bioluminescence imaging, 100 μL of D-luciferin (Yeason, Shanghai, China) was injected into each mouse and images were subsequently captured. The results were then analyzed using the Living Image software (Xenogen Corp.).

IHC staining

Three NPC and three adjacent normal paraffin-embedded tissues were collected from Guangdong Provincial People’s Hospital. All patients were pathologically confirmed. This study was approved by the Ethics Committee of Guangdong Provincial People’s Hospital. Tissue sections underwent deparaffinization using xylene followed by rehydration with a series of ethanol concentrations (100%, 95%, 85%, and 75%). Antigen retrieval was conducted before overnight incubation with the primary antibody at 4°C to inhibit endogenous peroxidase activity. Subsequently, the sections were incubated for 20 minutes with an HRP-conjugated secondary antibody at room temperature, and staining was carried out using diaminobenzidine (DAB) substrate (Dako). Following DAB treatment, the sections were counterstained with hematoxylin.

Statistical analysis

Statistical analyses within this research were executed using R and GraphPad Prism software. For comparisons between two groups, Student’s t-test and Wilcoxon test was used. When comparing more than two groups, one-way ANOVA and Kruskal-Wallis test was applied. Cox regression, log-rank tests were used to assess the prognostic significance of ELAVL1. Correlations between variables were assessed using Spearman’s correlation for non-parametric data and Pearson’s correlation for parametric data. The Benjamini–Hochberg method was applied to adjust the false discovery rate (FDR), with an adjusted P-value < 0.05 considered significant.

Expression patterns of ELAVL1 in pan-cancer and its localization

The human protein atlas (HPA) and FANTOM5 datasets revealed that ELAVL1 expression was most prominent within the thymus, testis, and lymph node tissues (Figure 2A and Supplementary Figure S1A). We also gathered data regarding ELAVL1 expression across various pan-cancer cell lines and found that it is highly expressed in adrenocortical cancer, leukemia, and lymphoma (Supplementary Figure S1B). Subsequently, we conducted a comparative analysis of ELAVL1 mRNA expression levels between tumor and normal tissues using data compiled from the TCGA and GTEx datasets, and we observed significantly increased expression of ELAVL1 in most cancer types, except for LAML, PCPG, and TGCT (Figure 2B and Supplementary Figure S1C). Similarly, ELAVL1 exhibited elevated mRNA expression levels in tumors compared to paired normal tissues, based on the TCGA pan-cancer dataset (Figure 2C and Supplementary Figure S1D). To verify its abnormal expression at the protein level, we collected pan-cancer proteomics datasets, and we noted that in most tumor types, there was a general upregulation of ELAVL1 protein levels (Figure 2D).

Figure 2

Expression patterns of ELAVL1 in pan-cancer and its localization. (A) The expression levels of ELAVL1 based on HPA database. (B) Violin plots showing the different mRNA expression levels of ELAVL1 in TCGA-GTEx pan-cancer datasets. (C) Paired differential analysis showing the mRNA expression levels between tumor and normal tissues in TCGA pan-cancer dataset. (D) Boxplots showing the different protein expression levels of ELAVL1 in proteomics pan-cancer datasets. (E) 3D structure of ELAVL1 obtained from PDB database. (F) Subcellular locations of ELAVL1 from the UniProt database. (G) Bar plots showing the AUC values evaluating the diagnostic efficacy of ELAVL1 expression in tumor and normal tissues. (H) Heatmap showing the mRNA expression level of ELAVL1 in pan-cancer datasets at the spatial transcriptome level. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001, and ^****P < 0.0001. HPA, human protein atlas; 3D, three-dimensional; PDB, Protein Data Bank.

The three-dimensional (3D) configuration of ELAVL1 was retrieved from the Protein Data Bank (PDB) database (Figure 2E). We then explored the cellular location of ELAVL1. Using the UniProt database, we found that ELAVL1 could be detected in the nucleoplasm, nucleoli, and cytosol (Figure 2F). The emergence of multi-omics technologies enables a more accurate exploration of the expression patterns of genes.^[48] To identify the particular cell types expressing ELAVL1, we examined its expression levels using the TISCH database-which contains pan-caner single-cell RNA sequencing (scRNA-seq) datasets-and the results showed that ELAVL1 was mainly expressed in malignant cells, proliferating T cells, and fibroblasts (Supplementary Figure S1E). To verify our finding, we collected pan-cancer spatial transcriptomics (ST) datasets. As expected, ELAVL1 was predominantly localized within the regions occupied by tumor cells (Figure 2H). In summary, our findings indicate that ELAVL1 is overexpressed in various cancers, which may underscore its potential as a clinical diagnostic biomarker.

ELAVL1’s provision of both diagnostic and prognostic value in different cancers

Receiver operating characteristic (ROC) curve analyses confirmed that ELAVL1 demonstrated high accuracy in distinguishing tumor and normal tissues for most cancer types such as BLCA, BRCA, and CESC (Figure 2G). Moreover, we performed univariate logistic regression, and the results showed that the odds ratio (OR) values were over 1 in most cancer types, suggesting a positive correlation between elevated ELAVL1 expression and tumor status and a possible association between increased ELAVL1 expression and an elevated risk of tumor development (Figure 3A).

To evaluate the prognostic utility of ELAVL1 across pan-cancer, we conducted Cox regression and log-rank tests for survival analyses encompassing overall survival (OS), disease-specific survival (DSS), progression-free interval (PFI), and disease-free interval (DFI) (Figure 3B). OS analysis showed that high ELAVL1 expression was associated with shorter OS in KIRC, KIRP, LGG, MESO, and SKCM, whereas it was a protective factor in CESC, KICH, LAML, OV, THCA, and THYM (Figure 3C). ELAVL1 was also a risk factor significantly associated with worse DFS in KIRC, KIRP, LGG, LUAD, MESO, and SKCM, while it was a protective factor in HNSC, KICH, and OV (Figure 3D). For PFI, ELAVL1 was a risk factor for KIRC, KIRP, and LGG, while it was protective in GBM and KICH (Figure 3E). With respect to DFI, ELAVL1 was a risk factor in BRCA and PRAD (Figure 3F). Overall, these findings indicate that ELAVL1 provides both diagnostic and prognostic value in different cancer types.

Analyses of ELAVL1 genomic alterations

Genomic approaches are a powerful means of cancer analysis.^{[49, 50]} Single nucleotide variation (SNV) mainly refers to the change in a DNA sequence caused by the alteration of a single nucleotide at the genomic level, which plays a significant role in tumor progression.^[51,52] Abnormal copy number variation (CNV) is also widely recognized as a fundamental molecular mechanism in tumor progression.^{[53, 54]} ELAVL1 amplification was chiefly detected in SARC, whereas deep deletions were frequently identified in UCEC, DLBC, SKCM, and SNVs were notably frequent in UCEC, DLBC, and STAD (Figure 4A–4D). We then explored CNV patterns in pan-cancer (including homozygous amplification, heterozygous amplification, homozygous deletion, and heterozygous deletion) (Figure 4D). Our analysis revealed that the CNV profiles of ELAVL1 predominantly exhibited heterozygosity (Figure 4E). Moreover, we analyzed the correlation of CNV with mRNA expression and found that most cancer types showed a positive correlation (Figure 4F).

Figure 3

ELAVL1 offers diagnostic and prognostic utility in different cancers. (A) Univariate logistic regression analysis of ELAVL1 mRNA expression and tumor status. (B) Heatmap showing relationship between ELAVL1 mRNA expression and different survival outcomes in pan-cancer. (C–F) Forest plots were used for pan-cancer analyses of ELAVL1 and OS (C), DFS (D), PFI (E), and DFI (F). OS, overall survival; DSS, disease-specific survival; PFI, progression-free interval; DFI, disease-free interval.

Figure 4

Analyses of ELAVL1 genomic alterations. (A) Stacked bar plots showing mutation frequencies of ELAVL1 in pan-cancer. (B) Heatmap showing mutation of ELAVL1 and several classical carcinogenic signaling pathways in pan-cancer. (C) Oncoplot of the mutation distribution of ELAVL1 in pan-cancer. (D) Percentage pie charts showing CNV profile of ELAVL1 in each tumor type. (E) The heterozygous and homozygous CNV profile of ELAVL1 in each tumor type, including the percentage of amplification and deletion. (F) Bubble plots showing the correlation of CNV with ELAVL1 mRNA expression. FDR, false discovery rate. CNV, copy number variation.

Modification patterns of ELAVL1 in pan-cancer

Epigenetic modifications and post-transcriptional modifications (PTMs) are crucial in influencing the onset and progression of cancers. ^{[55, 56, 57, 58]} By utilizing the UALCAN database, we first examined the promoter methylation levels of tumor and normal tissues, and the results showed that the ELAVL1 promoter methylation levels were upregulated in COAD, HNSC, KIRC, LUSC, and PAAD, whereas they were downregulated in PCPG (Supplementary Figure S2A). In addition, we observed a negative correlation between ELAVL1 mRNA expression levels and methylation levels (Supplementary Figure S2B).

We then explored the PTMs of ELAVL1. Figure 5A shows the common PTM sites of ELAVL1 obtained from the PhosphoSite database, indicating that S197, S200, and S202 were the most frequent phosphorylation sites. These findings were subsequently verified using the UACLAN database. We found that the ELAVL1 phosphorylation levels at these specific sites were elevated in tumor tissues (Figure 5B). Collectively, these results indicate the potential involvement of ELAVL1 in epigenetic modifications and PTMs.

Figure 5

Modification and alternative splicing patterns of ELAVL1 in pan-cancer. (A) An overview of Modification sites of ELAVL1 obtained from PhosphoSite database. (B) Box plots showing the different phosphorylated ELAVL1 levels between tumor and normal tissues using UALCAN database. (C) Reads-in, reads-out, and PSI values of ELAVL1_alt_3prime_134894 among pan-cancer, adjacent samples, and healthy tissues using OncoSplicing database. (D) PSI differences when comparing tumors and corresponding healthy or adjacent tissues and the association between ELAVL1_alt_3prime_134894 events and prognosis using OncoSplicing database. ^*P < 0.05, ^**P < 0.01, and ^****P < 0.0001.

Link of alternative ELAVL1 splicing to survival outcomes

Specific alterations in gene splicing within tumor cells might facilitate tumor progression, and discovering these alternative splicing (AS) events is important for enhancing survival outcomes. ^{[59, 60, 61]} The detection of clinically relevant ELAVL1 AS events was achieved using the OncoSplicing server. Figure 5C illustrates the percent spliced in (PSI) of the ELAVL1_alt_3prime event, designated 134894, and Figure 5D presents the PSI differences between tumor tissues and normal or adjacent tissues, as well as their relationships with prognosis. Although minimal differences in expression levels were observed, we found that ELAVL1 AS events might be correlated with the clinical outcomes of patients with ACC, KICH, LGG, UCEC, etc.

Involvement of ELAVL1 in multiple oncogenic pathways and its correlation with immune activity

A functional enrichment analysis was performed to clarify the role of ELAVL1 in tumor cells. We collected hallmark gene sets and performed GSEA. The result indicated that ELAVL1 was positively associated with several proliferation biological processes (e.g., G2M checkpoint, mitotic spindle, and DNA repair), and negatively related to classic inhibitory oncogenic pathways (e.g., p53 pathways and TNF-α pathways) and immune activity (e.g., interferon responses and complement) (Figure 6A). We then collected CRISPR screening data for ELAVL1 among pan-cancer cell lines and found a negative score for ELAVL1 across a substantial number of cell lines, indicating inhibition of cell proliferation and/or cell death following ELAVL1 knockout (Figure 6B).

Subsequently, we assessed the correlation between ELAVL1 expression levels and both genomic status and immune cell infiltration based on the immunogenicity and DNA damage scores. The results were mainly consistent with our previous findings, according to the quartile values of ELAVL1 expression levels (Figure 6C). In addition, we investigated the potential distinction among immune subtypes based on ELAVL1 expression levels and found that a high ELAVL1 expression level had a greater proportion of C1 subtype (increasing angiogenic gene expression and a higher proliferation percentage) but a lower proportion of C3 subtype (fewer aneuploidies and fewer CNV changes) (Figure 6D). The tumor microenvironment (TME) is vital in regulating tumor progression and modulating the response to standard-of-care therapies. ^{[62, 63, 64, 65]} Therefore, we additionally assessed the correlation between ELAVL expression levels and TME components and observed an immunosuppressive tendency (Supplementary Figure S3). Moreover, inverse negative correlations between ELAVL and immune-related genes existed in most cancer types, except THCA (Supplementary Figure S4). Pan-cancer immunotherapy cohorts were also summarized and the results indicated that ELAVL might felicitously predict immunotherapy response (Figure 6E).

ELAVL1 as a promising biomarker of NPC progression

Since no studies have examined ELAVL1 expression or its effect on NPC prognosis, we then focused on whether ELAVL1 could be a promising biomarker of NPC progression. We first collected four GEO datasets (GSE12452, GSE53819, GSE61218, and GSE103611) containing NPC and normal tissues. Compared to normal tissues, NPC tissues exhibited significantly increased ELAVL1 mRNA expression, and ELAVL1 was elevated in metastatic NPC tissues (Figure 7A). The AUC values for these datasets were 0.884, 0.778, 0.950, and 0.726, respectively, demonstrating the diagnostic efficacies of elevated ELAVL1 expression in NPC (Figure S6A). Besides, we found that elevated ELAVL1 mRNA expression levels were correlated with higher pathological N, stage, grade, and more specific mutations in the TCGAHNSC dataset (Figure S6B and S6C). We also observed upregulation of ELAVL1 mRNA expression levels in the NPC cell line with radiotherapy resistance (Figure 7B).

We then stratified the patients derived from the GSE102349 dataset into two groups characterized by the optimal cutoff value. We found that higher ELAVL1 expression levels were linked to unfavorable progression-free survival in NPC patients (Figure 7C). Taken together, our results suggest that ELAVL1 might be an effective biomarker of NPC progression.

ELAVL1 knockdown inhibited NPC progression in vitro

Subsequently, we performed western blotting to verify the upregulation of ELAVL1 protein levels in NPC cell lines compared to the normal nasopharyngeal cell line N₂-Tert. We then selected SUNE-1 and HONE-1 cell lines to determine whether ELAVL1 knockdown could inhibit NPC progression (Figure 7D). Small interfering RNAs (siRNAs) targeting ELAVL1 mRNA were designed and validated in NPC cell lines (Figure 7E). CCK-8 assays confirmed that ELAVL1 suppression led to a significantly decreased viability of NPC cell lines (Figure 7F). Additionally, the results of the Transwell and wound healing assays confirmed that ELAVL1 suppression markedly hindered NPC cell line migration (Figure 7G, 7H, and 7J). We subsequently performed colony formation assays and found that ELAVL1 knockdown significantly led to a reduction in NPC cell proliferation and stemness (Figure 7I and 7J). Moreover, we performed immunohistochemical (IHC) staining, and the results demonstrated that ELAVL1 was highly expressed in NPC tissues compared to adjacent normal tissues (Figure 8A). These in vitro findings underscore that ELAVL1 knockdown could impede the proliferation and migration of NPC cell lines.

ELAVL1 knockdown suppressed NPC progression in vivo

We established popliteal lymph node metastatic and lung metastatic models using the HONE-1 NPC cell line to further substantiate the role of ELAVL1 in vivo. Short hairpin RNA (shRNA) targeting ELAVL1 was applied and verified through western blotting (Figure 8B). ELAVL1-KD and control NPC cells were administered via injection into the footpads or tail veins of BALB/C nude mice (Figure 8C). For the popliteal lymph node metastatic model, we performed bioluminescence imaging and HE staining and observed that ELAVL1 knockdown inhibited the metastatic spread of NPC cells from the primary tumor located in the footpad to the popliteal space, resulting in tumor-draining lymph nodes with reduced sizes (Figure 8D, 8E, 8G, and 8H). Meanwhile, the tumor volumes at the primary site exhibited a significant reduction in the ELAVL1-KD group (Figure 8F–8H). For the lung metastatic model, the intravenous injection of ELAVL1-KD NPC cells exhibited a reduced lung metastatic burden compared to the control group, as evidenced by a decrease in overt metastatic lesions (Figure 8I–8L). Moreover, in these metastatic models, ELAVL1 knockdown inhibited the infiltration of NPC cells into normal tissues. Collectively, our data suggest that suppressing ELAVL1 expression might limit the migratory behavior of NPC cells and hinder metastasis in vivo.