Differential gene expression analysis reveals signaling pathways of proliferation, angiogenesis and metastasis is associated with interferon suppression in advanced hypopharyngeal cancer cells

Introduction

Hypopharyngeal cancer (HPC) is a highly aggressive subtype of head and neck squamous cell carcinoma (HNSCC) developing from the epithelial lining of the hypopharynx. In the United States, HPC accounts for 3–5% of all head and neck malignancies and is often diagnosed at an advanced stage due to its asymptomatic early progression and the anatomical complexity of the hypopharynx (1,2). Risk factors, including tobacco use, excessive alcohol consumption and human papillomavirus infection, contribute markedly to its pathogenesis (3,4). HPC is clinically classified according to the TNM staging system; early-stage disease (stage I and II) is often localized, while in its advanced stages (stage III and IV), regional lymph node involvement and distant metastasis are observed, leading to poor prognosis (5). Due to the aggressive nature of late-stage HPC, treatment strategies involve a combination of surgery, radiotherapy and chemotherapy; however, five-year survival rates have been reported to be 15–45%, due to high recurrence and resistance to therapy (6).

The progression of HPC, at the molecular level, is driven by a complex network of signaling pathways that regulate tumor growth, angiogenesis and metastasis. A hallmark of HPC is aberrant activation of the EGFR pathway, promoting uncontrolled cell proliferation through downstream effectors such as Ras-Raf-MEK-extracellular signal-regulated kinase and phosphoinositide 3 kinase-mammalian target of rapamycin, which promote tumor cell survival, proliferation and resistance to apoptosis (7). The epithelial-mesenchymal transition (EMT) process, mediated by transforming growth factor-β, Wnt/β-catenin and Notch signaling, is pivotal in enhancing tumor cell motility and invasiveness (8). Therefore, angiogenesis is a key factor supporting metastasis, driven by hypoxia-inducible factor-1α (HIF-1α) and vascular endothelial growth factor for neovascularization, ensuring an adequate supply of oxygen and nutrients to proliferating tumor cells (9). Increased angiogenic activity in advanced HPC is associated with increased metastatic potential and worse clinical outcomes (10). Recent studies have highlighted the role of HIF-1α in promoting metabolic reprogramming and angiogenesis, contributing to tumor progression and resistance to therapy (11,12).

Suppression of interferon (IFN) signaling has been identified as a key mechanism of immune evasion in HPC, reducing the susceptibility of the tumor to immune surveillance. IFN-stimulated genes (ISGs) are required for the enhancement of immune recognition and antiviral responses, while they have been reported to be frequently downregulated, impairing antitumor immunity (13,14). The loss of IFN signaling not only reduces immune-mediated tumor clearance but also promotes resistance to immunotherapy.

The present study aimed to identify the differentially expressed genes (DEGs) between parental FaDu HPC cells and advanced FaDu (FaDuex) cells, which were isolated from a late-stage xenograft tumor, were analyzed using the RNA-sequencing (RNA-seq) technique.

Materials and methods

Cell lines

The human HPC FaDu cells (cat no. HTB-43TM; American Type Culture Collection) and FaDuex cells were cultured in RPMI-1640 medium (Thermo Fisher Scientific, Inc.) containing 10% FBS (Thermo Fisher Scientific, Inc.), 1% penicillin/streptomycin and 1% L-glutamine (both MilliporeSigma; Merck KGaA). The FaDu cells were originally derived from a punch biopsy of a hypopharyngeal tumor obtained from a 56-year-old White male patient (15). FaDuex cells were previously established and considered as late-stage or advanced HPC cells, isolated from FaDu cells-derived xenograft tumors that reached 2,000 mm3 at 40 days (10). In brief, when the tumor size reached 600 mm3 exhibiting a necrotic region, the tumor was excised and cultured in RPMI-1640 medium supplemented with high-dose antibiotics, allowing viable cells to migrate from the tumor mass. These cells have been demonstrated to be more tumorigenic, exhibiting an increased rate of proliferation, angiogenesis and invasive capacity (10). The biological characteristics of FaDuex cells closely reflect the clinical features of human HPC diagnosed at a late stage (16).

RNA extraction and RNA-seq analysis

RNA was extracted from cells using the GENEZol™ TriRNA Pure Kit (cat. no. GZX100; Geneaid Biotech), according to the manufacturer's instructions. Up to 5×106 cells were harvested by centrifugation at 300 × g for 5 min at room temperature and the culture medium was completely removed. The cell pellet was thoroughly mixed in 700 µl of GENEzol™ reagent by pipetting and incubated for 5 min at room temperature. Ethanol was added in a 1:1 ratio to the lysate and the mixture was processed through a spin column-based protocol. RNA was eluted using RNase-free water and stored at −80°C until further use. For each condition (FaDu and FaDuex), there were two biological replicates and RNA was independently extracted from separate cultures, and reproducibility was ascertained. RNA library preparation and next-generation sequencing were performed by Biotools Co., Ltd. using their RNA-Seq (Q)20M package (service code: B-IRQTNT-E20P). Libraries were constructed using poly(A) mRNA enrichment and sequenced with 150 bp paired-end reads on the Illumina platform. All samples were processed and sequenced in a single batch.

The raw sequencing reads were processed for quality control using Trimmomatic (version 0.38) by Biotools Co., Ltd. and assessed using FastQC and MultiQC (17,18). Clean reads were aligned using HISAT2 to the human reference genome (GRCh38) and read counts for each gene were quantified using ‘featureCounts’. DEGs were analyzed by comparing the transcripts per million (TPM) values of genes between FaDuex and FaDu samples using ‘DESeq2’ to calculate fold-changes and statistical significance of the outcomes (19). Gene set enrichment analysis (GSEA) in ‘clusterProfiler’ (version 4.6.2, Bioconductor release 3.16) (20,21) was performed to identify significant pathways in the Hallmark and C2 gene sets from the Human Molecular Signatures Database (MSigDB; h.all.v2023.2.Hs and c2.all.v2023.2.Hs) (22–24), as well as other pathways of interest, including 23 downstream NFE2 like BZIP transcription factor 2 gene sets [Transcription_Factor_PPIs (25), Rummagene_transcription_factors (26), TRRUST_Transcription_Factors (27) and TRANSFAC_and_JASPAR_PWMs (28)] (26,28,29) and 20 cancer stem cell-related gene sets (30) collected from Enrichr (25,31,32). Only pathways with an adjusted P<0.05, as determined by Benjamini-Hochberg correction, were considered.

Complementary DNA (cDNA) synthesis and quantitative polymerase chain reaction (qPCR)

The cDNA was synthesized using the ToolsQuant II Fast RT kit (cat. no. KRT-BA-06; TOOLS Biotech) following the manufacturer's instructions. Briefly, after quantifying total RNA using a NanoDrop ND-1000 Spectrophotometer (Thermo Fisher Scientific, Inc.), 300 ng of RNA was used as the template for reverse transcription. Genomic DNA contamination was removed by incubating with gDNA Eraser at 42°C for 3 min, followed by reverse transcription at 42°C for 15 min. The resulting cDNA was diluted with 80 µl nuclease-free water and utilized directly for subsequent qPCR analysis employing the StepOnePlus™ Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). Primers used for qPCR are summarized in Table I. Each reaction was carried out in a final volume of 20 µl, comprising 6 µl nuclease-free water, 2 µl primer mix (6 µM forward primer and 6 µM reverse primer), 2 µl cDNA template and 10 µl, EnTurbo™ SYBR Green PCR SuperMix (cat. no. EQ013; High ROX Premixed; ELK Biotechnology). The thermal cycling program was as follows: Initial denaturation at 95°C for 20 sec, followed by 40 cycles of denaturation at 95°C for 3 sec and annealing/extension at 60°C for 30 sec, during which fluorescence signals were recorded. At the end of the amplification reaction, the melt curve was analyzed to confirm product specificity. The expression of each target gene, along with the internal control gene (S26), was analyzed in four technical replicates. Each datum represented a mean of four repeats and was compared using unpaired t-tests. Relative expression levels were estimated based on amplification curves (33).

Table I. Sequences of primers for each gene analyzed by quantitative PCR.

Table I.

Sequences of primers for each gene analyzed by quantitative PCR.

Gene	Primer sequence (5′-3′)
KRT13	F: GATGCTGAGGAATGGTTCCACG
	R: AGCTCCGTGATCTCTGTCTTGC
IFI44L	F: TGCACTGAGGCAGATGCTGCG
	R: TCATTGCGGCACACCAGTACAG
IFI6	F: TGATGAGCTGGTCTGCGATCCT
	R: GTAGCCCATCAGGGCACCAATA
STC1	F: GCAGGAAGAGTGCTACAGCAAG
	R: CATTCCAGCAGGCTTCGGACAA
MUC16	F: GATGTCAAGCCAGGCAGCACAA
	R: GAGAGTGGTAGACATTTCTGGGC
S26	F: CCGTGCCTCCAAGATGACAAAG
	R: ACTCAGCTCCTTACATGGGCTT

[i] KRT13, keratin 13; IFI44L, interferon-induced protein 44-like; IFI6, interferon α-inducible protein 6; STC1, stanniocalcin-1; MUC16, mucin-16; F, forward; R, reverse.

Results and Discussion

Analysis of DEGs between parental and advanced HPC cells isolated from late-stage xenograft tumors

The present study first analyzed the DEGs between FaDu and FaDuex cells (advanced FaDu cells) to identify transcriptional differences associated with their phenotypic variation. The heatmap in Fig. 1A presents distinct clustering patterns; genes upregulated in FaDuex cells are displayed in red and the downregulated genes are presented in green, indicating significant transcriptional reprogramming (Fig. 1A). The differences are further illustrated in the volcano plot (Fig. 1B). A total of 86 genes were significantly upregulated, while 166 genes were significantly downregulated in FaDuex cells relative to their expression in FaDu cells (adjusted P<0.05). These findings suggest that the expression levels of genes in FaDuex cells undergo significant changes, potentially explaining their distinct biological properties.

Figure 1. Comparison of DEGs between FaDuex and FaDu cells. (A) Heatmap of 252 DEGs between FaDuex and FaDu cells, with log2 fold-change ≥1.5 or ≤-1.5 and adjusted P<0.05; (B) the volcano plot displays upregulated and downregulated genes, with the x-axis representing the log2 fold-change and the y-axis representing the -log10 adjusted P-value with log2 fold-change ≥1.5 or ≤-1.5 and -log10(adjusted P-value) >-log10(0.05). DEGs, differentially expressed genes. TPM, transcripts per million; FaDuex, advanced FaDu.

Comparison of specific signaling pathways activated and suppressed in advanced HPC cells

GSEA was performed next using the Hallmark and C2 gene sets to investigate the biological processes associated with the DEGs between FaDuex and FaDu cells. Significantly enriched pathways were subsequently categorized into activated and suppressed groups. In the Hallmark gene set, only suppressed enriched pathways were significant and the most suppressed pathways included those involved in ‘IFNα response’, downregulated ‘KRAS signaling’, ‘IFNγ response’, ‘myogenesis’ and ‘heme metabolism’ (Fig. 2A). In the C2 gene set with significant changes, four enriched pathways (upregulated ‘Wotton RUNX targets’, ‘AMIT EGF response’ in MCF10A and HeLa cells and ‘Pedersen metastasis by ERBB2 isoform’) were activated, while 10 pathways were suppressed, including those involved in ‘BROWNE IFN response’ and ‘Pedersen IFN-α response’ (Fig. 2B). Since the IFN pathways are key for antiviral defense and tumor immune surveillance (34,35), their suppression in FaDuex cells possibly indicates an immune evasion phenomenon, potentially facilitating tumor progression and resistance to immune-mediated clearance. These findings highlighted a significant suppression of IFN signaling in FaDuex cells, potentially contributing to their altered immune landscape and tumorigenic potential.

Figure 2. Gene set enrichment analysis between FaDuex and FaDu cells, using the human molecular signatures database. The image displays (A) top 10 significantly suppressed signaling pathways in FaDuex cells, identified using the hallmark gene set. (B) top signaling pathways activated and suppressed in FaDuex cells using the C2 gene set. The number of genes involved is represented by the size of each dot and the color indicates the adjusted P-value. FaDuex, advanced FaDu.

Gene regulations of advanced HPC cells involved in different cancerous associated signaling pathways

The changes in the expression levels of genes associated with several pathways involved in metastasis, angiogenesis, EGF signaling and EMT in FaDuex compared with FaDu cells were further examined. Assessment of differential expression of genes revealed that several metastasis-related genes were altered, with keratin 13 (KRT13), RIBC2 and serpin family B member 13 being the most significantly downregulated metastasis-inhibiting genes, and stanniocalcin-1 (STC1), high mobility group nucleosome binding domain 5 and down syndrome cell adhesion molecule being the most significantly upregulated metastasis-promoting genes, suggesting an increase in metastatic potential in FaDuex cells (Fig. 3A). Regarding the angiogenesis pathway, only interferon α-inducible protein 6 (IFI6), insulin-like growth factor binding protein 7 and intercellular adhesion molecule 1 (ICAM1) were significantly downregulated angiogenesis-inhibiting genes in FaDuex cells, indicating a possible promotion of neovascularization (Fig. 3B). Furthermore, a significant upregulation of STC1, doublecortin-like kinase 1, transglutaminase 2, cytochrome P450 family 27 subfamily B member 1 and forkhead box D1 genes was noted, indicating the activation of the EGF-related pathway, which may contribute to enhanced proliferative and survival signaling in FaDuex cells (Fig. 3C). Furthermore, the genes involved in the inhibition of the EMT pathway, including KRT13, interleukin 1 receptor antagonist, aquaporin-3, complement component 1s, tripartite motif-containing protein 29 and transforming growth factor-β receptor 3, were downregulated in FaDuex cells, suggesting a shift toward EMT phenotypes (Fig. 3D). These findings indicated that FaDuex cells exhibit enhanced metastatic, angiogenic, proliferative and EMT-related potential, reflecting a more aggressive phenotype.

Figure 3. Heatmaps of Z-score normalized TPM values and bar plots of log2 fold-change of DEGs. The DEGs were associated with (A) metastasis, (B) angiogenesis, (C) proliferation and survival and (D) epithelial-mesenchymal transition-related genes. The color gradient indicates the Z-score of TPM. Bars were colored black for adjusted P≥0.05 and log2 fold-changes between −1.5 and 1.5. For bar plots with adjusted P<0.05, log2 fold-changes <-1.5 are marked in green, while those >1.5 are marked in red, indicating significant fold-changes. TPM, transcripts per million; DEGs, differentially expressed genes; FaDuex, advanced FaDu.

Downregulation of interferon signaling-associated genes in advanced HPC cells

The present study analyzed the expression of IFN-responsive genes using the C2 and Hallmark gene sets to evaluate the differences in IFN signaling between FaDu and FaDuex cells. Both analyses revealed a broad suppression of ISGs in FaDuex cells, indicating a significant suppression of IFN signaling. In the C2 gene set (Fig. 4A), key ISGs such as interferon-induced protein 44-like (IFI44L), IFI6 and tripartite motif-containing protein 22, all of which are key mediators of the IFN response and antiviral defense mechanisms, were markedly downregulated. The Hallmark gene set (Fig. 4B) further confirmed the suppression of IFN signaling, demonstrating a decrease in the expression of genes such as IFI44L, caspase-1 and ICAM1, which are involved in immune activation and inflammatory responses. In the Hallmark gene set, mucin-16 (MUC16) was the most significantly downregulated; however, its expression was only stimulated by the combined action of tumor necrosis factor-α and IFN-γ (36). The downregulation of these ISGs suggested that FaDuex cells may have reduced sensitivity to IFN-mediated immune responses, potentially altering their interaction with the immune microenvironment and affecting their ability to respond to external immune challenges.

Figure 4. Heatmaps of Z-score-normalized TPM values and bar plots of log2 fold-changes for differentially expressed genes associated with the interferon-responsive signaling pathway. (A) Analysis of interferon-stimulated genes using the C2 gene set and (B) the hallmark gene set. The color gradient indicates the Z-score of TPM. Bars are colored black for adjusted P-values ≥0.05 and log2 fold-changes between −1.5 and 1.5. For bar plots with adjusted P-values <0.05, log2 fold-changes <-1.5 are marked in green, while those >1.5 are marked in red, indicating significant fold-changes. TPM, transcripts per million; FaDuex, advanced FaDu.

Validation of gene expression using qPCR

The results of DEGs analysis were validated by performing qPCR by selecting genes, in FaDu cells and FaDuex cells, whose expression levels exhibited the most significant changes. According to the DEGs analysis, KRT13 and IFI6 were most significantly downregulated in metastasis/EMT and angiogenesis, respectively, whereas STC1 was the most significantly upregulated in proliferative and survival-related gene of FaDuex cells. In addition, IFI44L and MUC16 were most downregulated in the IFN-responsive genes of FaDuex cells according to the C2 and Hallmark gene sets. The qPCR analysis further demonstrated the downregulation of KRT13, IFI44L, IFI6 and MUC16 transcripts (Fig. 5A, B, D and E) and STC1 was upregulated in FaDuex cells (Fig. 5C), consistent with DEG analysis. These results suggested that RNA-seq determined gene expression is useful for the prediction of cancer-related signaling pathways.

Figure 5. Quantitative PCR analysis of genes expressed in FaDu cells and FaDuex cells. (A) KRT13; (B) IF16; (C) STC1; (D) IFI44L and (E) MUC16. Data are presented as mean ± SD from four independent experiments. **P<0.01 and ****P<0.0001. KRT13, keratin 13; IFI44L, interferon-induced protein 44-like; IFI6, interferon α-inducible protein 6; STC1, stanniocalcin-1; MUC16, mucin-16; FaDuex, advanced FaDu.

The limitations of the present study include a small sample size; in addition, the results were analyzed from human HPC cell lines rather than clinical samples. As mentioned, HPC is a rare subtype of HNSCC (6). In Taiwan, the age-standardized incidence rate for HPC was increased from 1980 and reached 6.46 per 100,000 in 2019 (37). A pilot study may be initiated in the future to determine notable changes of gene expression using a widely used HPC cell line and its derivatives from xenograft tumor. FaDu cells are isolated from the primary tumor of a White male patient with HPC in 1968 (15). Recently, another HPC cell type (CZH1) isolated from a Chinese patient with HPC has been reported, exhibiting a greater capacity for invasion and increased radiosensitivity compared with FaDu cells (38). The present study would be key in exploring the gene expression differences of patients with HPC from different ethnicities. FaDuex cells are isolated from xenograft tumors that have reached 2,000 mm3 and the present study defined them as late-stage tumors because of increased metastatic, tumorigenic, angiogenic and chemoradiotherapy-resistant capacity (10). However, it is difficult to explain the TNM stage of clinical HPC through late-stage FaDuex cells, as they exhibit similar characteristics as advanced human HPC diagnosed in clinics. In addition, the translational limitation may still exist when using FaDuex cells to represent late-stage HPC in clinical settings. Another limitation of the present study was the use of parental FaDu cells for comparison with tumor-derived FaDuex cells rather than the early stage of FaDu tumors (for example, tumor size at 100 mm3). It may be key to compare potent DEGs in FaDu-derived xenograft tumors at the early (small size) stage and the late (large size) stage in the future.

To conclude, the present study demonstrated that FaDuex cells, the arbitrarily defined late-stage HPC cells isolated from xenograft tumors, exhibited several genes - including KRT13, IFI6, STC1, MUC16 and IFI44L - that were differentially expressed and were associated with malignant characteristics, including increased metastasis, angiogenesis, proliferation and decreased inflammatory responses. The present study also observed that RNA-seq data-based DEGs analysis could be validated using qPCR. Although the use of human cell lines may encounter several of the aforementioned limitations, the variation in samples would be low and easy to interpret. The present study results suggest that several genes are associated with the advanced malignancy of HPC cells, which makes this information potentially key for diagnostic and therapeutic considerations in the future.

Acknowledgements

Not applicable.

Funding

The present study was supported by a grant of Taipei City Hospital, RenAi Branch (grant no. TPCH-112-12) and that of National Science and Technology Council (grant no. NSTC 114-2314-B-A49-065-MY3).

Availability of data and materials

The data generated in the present study may be found in the Sequence Read Archive (SRA) of National Library of Medicine under accession no. PRJNA1282043 or at the following URL: https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA1282043.

Authors' contributions

CWK conducted the cell culture experiments, RNA preparation, RNA-sequencing arrangement and data analysis. YPY conducted the differential gene expression and gene set enrichment analyses. TWC, JDL and YJL confirmed the authenticity of all the raw data. TWC interpreted the RNA sequencing data, and resolved questions related to the accuracy of this work. MYL isolated and identified different types of FaDu cells. JDL and YJL contributed to the conception and the design of the present study. YJL drafted the manuscript. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Use of artificial intelligence tools

During the preparation of this work, AI tools were used to improve the readability and language of the manuscript or to generate images, and subsequently, the authors revised and edited the content produced by the AI tools as necessary, taking full responsibility for the ultimate content of the present manuscript.

References