Genomic and Single-Cell Analyses Characterize Patient-Derived Tumor Organoids to Enable Personalized Therapy for Head and Neck Squamous Cell Carcinoma

Abstract

Head and neck squamous cell carcinoma (HNSCC) remains a significant health burden because of tumor heterogeneity and treatment resistance, emphasizing the need for improved biological understanding and tailored therapies. In this study, we enrolled 31 patients with HNSCC for the establishment of patient-derived tumor organoids (PDO), which faithfully maintained the genomic features and histopathologic traits of the primary tumors. Long-term culture preserved key characteristics, affirming PDOs as robust representative models. PDOs demonstrated predictive capability for cisplatin treatment responses, with ex vivo drug sensitivity correlating with patient outcomes. Bulk and single-cell RNA sequencing unveiled molecular subtypes and intratumor transcriptional heterogeneity (ITH) in PDOs, paralleling patient tumors. Notably, a hybrid epithelial-mesenchymal transition-like ITH program was associated with cisplatin resistance and poor patient survival. Functional analyses identified amphiregulin as a potential regulator of the hybrid epithelial-mesenchymal state. Moreover, amphiregulin contributed to cisplatin resistance via EGFR pathway activation, corroborated by clinical samples. In summary, HNSCC PDOs serve as reliable and versatile models, offer predictive insights into ITH programs and treatment responses, and uncover potential therapeutic targets for personalized medicine.

Significance: Profiling of patient-derived organoids uncovers intertumoral heterogeneity and a hybrid epithelial-mesenchymal transition program conferring cisplatin resistance and highlights amphiregulin as a regulator of cellular plasticity and potential therapeutic target for HNSCC treatment.

Figure 1.. HNSCC PDOs maintain histopathological characteristics of the original tumors.

(A) Workflow of the study design (created in BioRender. Nam, C. (2024) https://BioRender.com/v421oy3). (B) Oncoplots (left) of driver mutations and Venn diagrams (right) of total nonsynonymous mutations from two HNSCC PDO lines grown in different culture media. (C) H&E and IHC staining of indicated markers on selected PDOs and their matched primary tumors. Early, early-passaged PDOs; Late, late-passaged PDOs; MUT, mutant; WT, wildtype. (D) The expression of epithelial markers (E-cadherin and β-catenin) were evaluated using IF staining in selected PDO lines.

Figure 2.. Characterization of HNSCC PDO xenograft (PDOX).

(A) Images and (B) growth curves (mean±SD) of PDOX tumors. NT, no tumor formation. (C) H&E and IHC staining of matched PDOX, originating tumors and PDO samples. (D) IC50 values of each PDO line treated with cisplatin. PDOs were divided into sensitive and resistant groups based on the median IC50 (15 uM, indicated by the dashed line). (E) A Kaplan Meier survival plot of HNSCC patients, grouped according to the cisplatin sensitivity of their corresponding PDO lines.

Figure 3.. HNSCC PDOs preserve the genomic landscape of the parental tumors.

(A) Oncoplots of driver mutations and copy number variations (CNVs) comparing parental tumors (T) vs. early-passaged PDOs (E), or (B) between early- and late-passaged PDOs (L). TCGA-curated genomic pathways are annotated on the left. (C) Column plots summarizing the concordance rate of total nonsynonymous mutations between parental tumors vs. early-passaged PDOs, or (D) between early- and late-passaged PDOs. (E) Heatmap showing the unsupervised clustering of parental tumors and early-passaged PDOs. (F) CNV profiles of parental tumors and early-passaged PDOs. T: primary tumor, E: early-passaged PDOs, L: late-passaged PDOs.

Figure 4.. Bulk and single-cell transcriptomic analysis of HNSCC PDOs.

(A) Heatmap and (B) PCA plot showing molecular subtyping of 13 HNSCC PDOs based on bulk RNA-seq. (C) A UMAP showing 68,919 cells obtained from 8 PDO lines, colored by patient origin. (D) UMAP of the expression of selected marker genes. (E) Heatmap showing consistent molecular subtyping of PDOs using either scRNA-seq or bulk RNA-seq. (F) UMAP of subtype annotation of single cells. (G) Pearson correlation of subtype scores derived by scRNA-seq vs. bulk RNA-seq. (H) Pearson correlation of subtype scores derived by two independent methods based on scRNA-seq.

Figure 5.. scRNA-seq reveals the intratumor heterogeneity (ITH) of HNSCC PDOs.

(A) Heatmap showing unsupervised clustering of 50 NMF programs based on top program genes, with similarity measured by Jaccard indices. Eight clusters of consensus ITH programs are numbered and highlighted. (B) Dot plots showing the Jaccard index (y) and mean correlation (x) between ITH programs and published meta-programs (15). Dashed lines indicate a 99.9% confidence threshold determined using permutations of NMF programs. (C) A heatmap showing expression levels of program genes, with each program aligned to panel (A). (D) The most significantly correlated ITH programs from panel (B). (E-G) Heatmaps showing similarity (measured by Jaccard indices) between PDO ITH programs and published programs from primary HNSCC (E), ESCC (F) and CSCC (G) (14,33,34). *, FDR q < 0.05; **, FDR q < 0.01; ***, FDR q < 0.001; hypergeometric test. (H) PCA plots of PDO cells only (yellow, lower left), primary HNSCC tumor cells only (purple, lower middle) or combined (lower right). In the upper panel, cells are colored by the relative score for the hEMT minus EpiSen program identified from PDO samples. (I-J) similar PCA plots to panel H, except using program genesets identified from HNSCC primary tumors (I) or shared program genesets between PDO and primary tumors (J).

Figure 6.. The hEMT program is associated with PDO cisplatin resistance.

(A) UMAPs showing the distribution of each ITH program at the single-cell level. (B) Dot plots showing the program score and ratio of each ITH program in each PDO sample. (C-D) Box plots showing the program score (C) and cell ratio (D) of each ITH program comparing the cisplatin-sensitive group (n=4) with the resistant group (n=4). (E-F) Pearson correlation between program score (E) and cell ratio (F) of each ITH program vs. IC50 values of cisplatin. Each dot is one PDO sample. (G) A Kaplan meier survival plot of TCGA HNSCC patients (n=520), grouped according to the mRNA expression of 32 hEMT genes using the GSVA method. (H) multivariate survival analysis of indicated factors using the TCGA HNSCC cohort.

Figure 7.. AREG, a top hEMT program gene, regulates PDO responsiveness to cisplatin.

(A) A scatter plot showing the Pearson correlation R value (y) and p value (x) of each hEMT gene with IC50 values of cisplatin. Top 5 genes are highlighted. (B) Pearson correlation between AREG expression and IC50 values of cisplatin. Each dot is one PDO sample. (C) A column plot showing the expression fold change of seven EGFR ligands between HNSCC tumors with normal samples using TCGA bulk RNA-seq data. *, P < 0.05. (D) A Kaplan Meier survival plot of TCGA HNSCC patients (n=520), grouped according to the mRNA expression of AREG. (E) multivariate survival analysis of indicated factors using the TCGA HNSCC cohort. (F) A box plot showing the protein levels of AREG across various HPA human tumor samples. (G) Box plots showing the mRNA expression of AREG in indicated HNSCC tumor samples from the dataset GSE9349. (H) dose-response curves of scramble control and AREG-knockdown samples to cisplatin treatment in HNSCC cell lines (A253 and SNU-1066) or (I) PDOs (PDO-42 and PDO-49). (J) Images of PDOs upon cisplatin treatment. (K) dose-response curves and (L) images of PDO-51 to cisplatin treatment in the presence and absence of recombinant AREG protein. Scale bar=40 um

Figure 8.. AREG regulates the activity of hEMT program in HNSCC.

(A) Western blotting of indicated antibodies upon either knockdown of AREG with siRNAs or addition of recombinant AREG protein in SNU-1066 cells. (B-C) hockey stick plots showing the correlation of the mRNA expression of AREG with other genes, rank ordered by correlation coefficient. hEMT program genes are highlighted in red. Only genes with nonnegative correlation are shown. (B) is the TCGA dataset and (C) is CCLE+HPA combined dataset. (D-F) heatmaps showing changes in expression of hEMT program genes upon either AREG knockdown, or addition of recombinant AREG protein, in HNSCC cell lines or PDO samples. (G) A schematic graph of the experiment design (created in BioRender. Nam, C. (2025) https://BioRender.com/8hxao4r). (H) Growth curves and (I) images of xenograft tumors. (J) Representative IHC images of p-ERK and p-AKT on each group. Scale Bar=100 um.