To systematically identify GPCRs associated with pancreatic cancer prognosis, the present study integrated multiple genomic resources. GPCR-related genes were curated from WikiPathways (IDs: WP24, WP58, WP80, WP117, WP247, WP334, WP455, WP501; www.wikipathways.org), while prognostic genes and highly expressed transcripts in pancreatic adenocarcinoma were retrieved from UALCAN (https://ualcan.path.uab.edu) and GEPIA2 (http://gepia2.cancer-pku.cn), respectively (Table SI, Table SII, Table SIII). Transcriptomic profiles (TPM format), clinical information and somatic mutation data for 178 patients with pancreatic cancer were acquired from TCGA-PAAD cohort using the ‘TCGAbiolinks’ R package (version 2.36.0; bioconductor.org/). ‘TCGAvisualize_oncoprint’ function was employed to generate a heatmap illustrating gene mutation frequencies. For independent validation, RNA-seq data and clinical records from the PACA-AU cohort were retrieved from the ICGC database (https://dcc.icgc.org). Raw read counts were converted to TPM values using the ‘count2tpm’ function in the ‘IOBR’ package (version 0.99.9; github.com/). Additionally, two microarray datasets (GSE28735 and GSE183795) based on the GPL6244 platform were incorporated to assess the predictive performance of the model (14,15). Clinical metadata and RMA-normalized expression matrices were downloaded from the GEO data portal (https://www.ncbi.nlm.nih.gov/geo). To ensure data comparability, the present study applied rigorous harmonization procedures. All transcriptomic data were log2-transformed [log2(x+1)]. For batch effect correction, the present study used the ComBat algorithm as implemented in the sva R package (v3.56.0, Bioconductor). Specifically, microarray datasets (GSE28735 and GSE183795) were based on RMA-normalized expression matrices, while RNA-seq data from TCGA and ICGC were converted to TPM (Transcripts Per Million) format before log2 transformation and batch correction.

To evaluate the expression profiles of HRH1 in matched tumor and adjacent non-tumor tissue samples, bulk RNA-seq datasets were retrieved from the GEO and ArrayExpress (https://www.ebi.ac.uk/biostudies/arrayexpress) databases, including GSE183795 (85 pairs), GSE62452 (42 pairs), GSE15471 (35 pairs), GSE101448 (18 pairs), GSE16515 (15 pairs), GSE196009 (6 pairs) and E-MTAB-6690 (65 pairs) (14,16–21). Additionally, transcriptomic data from GSE26088 were analyzed to evaluate HRH1 expression across 19 pancreatic cancer cell lines and one normal pancreatic ductal epithelial cell line (HPDE). PL5, also called Panc 04.03, is a pancreatic cancer cell line (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM640500; http://www.cellosaurus.org/CVCL_1636). To investigate the regulatory relationship between YAP1 and HRH1 across diverse cancer types, transcriptomic datasets from the GEO were utilized, specifically: Human endothelial cells (GSE211726), breast cancer (GSE59232), neuroblastoma (GSE130401), liver cancer (GSE137915) and renal clear cell carcinoma (GSE146354) (22–25). For single-cell resolution analysis, the CRA001160 dataset containing both normal and pancreatic cancer tissues was employed to compare HRH1 expression across distinct cell populations. Single-cell RNA sequencing data were analyzed with Seurat (v5.3.0; Bioconductor) in R, employing analytical approaches consistent with established protocols (26). Gene expression levels were quantified based on mean unique molecular identifier (UMI) counts.

Tumor mutational burden (TMB) calculation. TMB was calculated using somatic mutation data derived from whole-exome sequencing. Specifically, the present study retrieved the harmonized masked somatic mutation data from TCGA database via the ‘TCGAbiolinks’ R package. To compute TMB, the present study first filtered the mutation annotation format file to retain only protein-altering somatic variants, including missense, non-sense, splice-site and frameshift indel mutations. Variants classified as germline, silent or located in non-coding regions were excluded. TMB was defined as the total number of qualifying somatic mutations divided by the effective coding region size. Given that TCGA exome capture kits cover ~38 Mb of the coding genome, the present study used this value as the denominator, consistent with prior TCGA-based TMB studies (27,28). TMB values are reported in units of mutations per megabase (mut/Mb), with TMB-high cases defined as those exhibiting ≥ 10 muts/Mb (29).

Two-sample Mendelian randomization analysis. ‘TwoSampleMR’ package (version 0.6.14) from the Comprehensive R Archive Network (CRAN, cran.r-project.org/) was used to investigate potential causal relationships between exposure to fexofenadine therapy (GWAS ID: ukb-b-10433; UK Biobank, http://www.ukbiobank.ac.uk/) and pancreatic cancer (GWAS ID: ebi-a-GCST90018893; GWAS Catalog, http://www.ebi.ac.uk/gwas/). Fexofenadine therapy served as the exposure variable. Instrumental variables (IVs) were initially selected using a significance threshold of P<4.5×10⁻⁵. To mitigate the effects of linkage disequilibrium, IV clumping was performed with an r² threshold of 0.001 and a clumping distance cut-off of 10,000 kilobases. Weak IVs, characterized by an F-statistic <10, were subsequently excluded. Given the established associations of obesity, smoking and alcohol consumption with pancreatic cancer risk, IVs potentially associated with these confounding factors were identified through association queries in the GeneAtlas database (geneatlas.roslin.ed.ac.uk/) and subsequently excluded (5). The inverse-variance weighted (IVW) method served as the primary Mendelian randomization (MR) analysis. Secondary MR analyses employed the following methods: MR-Egger regression, weighted median, weighted mode and simple mode. Pleiotropy was assessed using the MR-Egger intercept test. Sensitivity analyses included leave-one-out validation and Cochran's Q statistic to evaluate the heterogeneity among single-nucleotide polymorphisms (SNPs).

The non-cancerous control cell line hTERT-HPNE (Zhejiang Meisen Cell Technology Co., Ltd.) was cultured in complete DMEM at 37°C within a humidified atmosphere containing 5% CO₂. The present study obtained all cell culture media, FBS and PDAC cell lines from Procell Life Science & Technology Co., Ltd. Cell lines were cultured under standardized conditions: ASPC1 and BXPC3 in RPMI 1640, CFPAC1 in Iscove's Modified Dulbecco's Medium (IMDM), while PANC1, MIAPaCa-2, Capan-2 and SW1990 were propagated in DMEM.

To knock down target gene expression, transient transfection of siRNAs was carried out with Lipofectamine 3000^® (Invitrogen; Thermo Fisher Scientific, Inc.) as per the manufacturer's guidelines. Briefly, 20 µM siRNA duplex was mixed with 6 µl Lipofectamine 3000 in serum-free medium and incubated at room temperature for 15 min to form complexes. The mixture was then added dropwise to cells seeded in 6-well plates containing 2 ml complete growth medium. Cells were maintained at 37°C for 12 h, after which the transfection medium was replaced with fresh complete medium. Experiments were conducted 48 h post-transfection to allow sufficient knockdown of target genes. The following small interfering (si)RNA duplexes (Shanghai GenePharma Co., Ltd.) were employed: For HRH1 silencing, i) 5′-GGACAAGUGUGAGACAGACTT-3′ (sense) and 5′-GUCUGUCUCACACUUGUCCTT-3′ (antisense); ii) 5′-GCUCUGGUUCUAUGCCAAGAU-3′ (sense) and 5′-AUCUUGGCAUAGAACCAGAGC-3′ (antisense); for YAP1 knockdown, 5′-GUCAGAGAUACUUCUUAAATT-3′ (sense) and 5′-UUUAAGAAGUAUCUCUGACTT-3′ (antisense). A non-targeting control siRNA (5′-UUCUCCGAACGUGUCACGUTT-3′ (sense) and 5′-ACGUGACACGUUCGGAGAAUTT-3′ (antisense) served as a negative control.

Fexofenadine (cat. no. HY-B0801; MedChemExpress) was employed as a selective HRH1 antagonist, while histamine (cat. no. HY-B1204; MedChemExpress) served as the HRH1 agonist.

A seeding density of 2,000 cells/well was used in 96-well plates. Following incubation at 37°C for 2 h, the culture medium was replaced with Cell Counting Kit-8 (CCK-8; cat. no. GK10001; GlpBio Technology) working solution (medium: CCK-8=10: 1) and incubated for 3 h. Absorbance at 450 nm was measured at 24-h intervals for up to 96 h.

For cytotoxicity assessment, tumor cells were seeded at a density of 4,000 cells per well in 96-well plates and cultured at 37°C for 24 h prior to treatment with gemcitabine. (cat. no. GC16805; GlpBio Technology) at concentrations of 1×10⁻⁶, 1×10-³, 1×10-², 1×10-¹, 1, 10, 100, and 1,000 µM. Following 48 h of drug incubation at 37°C, optical density (OD) signals were measured. Cell viability was calculated with the formula: Cell viability=(OD of gemcitabine-exposed group-blank OD)/(OD of untreated control group-blank OD). The IC₅₀ was determined by non-linear regression analysis using GraphPad Prism (version 9.00; GraphPad Software, Dotmatics).

Total RNA was extracted from cells using the Seven RNAkey Reagent (cat. no. SM139-02; Sevenbio) and reverse transcribed into cDNA with the Sevenbio cDNA Synthesis kit (cat. no. SM135-02), following the manufacturer's protocol. Gene expression analysis was performed on a 7500 Fast Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) using SYBR Green qPCR MasterMix (Sevenbio). 36B4 (RPLP0) was selected as the internal reference gene for qPCR analysis because it exhibits stable and consistent expression in pancreatic cancer cell lines, which has been widely validated in pancreatic cancer gene expression studies (30–35). CTGF (also known as CCN2) and CYR61 (also known as CCN1) are key downstream effector genes in the YAP1 signaling pathway (36). All primer sequences provided below are presented from the 5′ to 3′direction. The primer sequences used for qPCR were as follows: 36B4 forward (F): GCAGCATCTACAACCCTGAAG; 36B4 reverse (R): CAC TGGCAACATTGCGGAC; CYR61 F: CCCGTTTTGGTAGATTCTGG; CYR61 R: GCTGGAATGCAACTTCGG; CTGF F: ACCGACTGGAAGACACGTTTG; CTGF R: CCAGGTCAGCTTCGCAAGG; ANKRD1 F: GTGTAGCACCAGATCCATCG; ANKRD1 R: CGGTGAGACTGAACCGCTAT; YAP1 F: CAGACAGTGGACTAAGCATGAG; YAP1 R: CAGGGTGCTTTGGTTGATAGT A; HRH1 F: GCTGGGCTACATCAACTCCAC; HRH1 R: CCCTTAGGAGCGAATATGCAGAA. The thermocycling conditions were 95°C for 30 sec (initial denaturation), followed by 40 cycles of 95°C for 10 sec (denaturation) and 60°C for 30 sec (annealing/extension). Fluorescence signals were collected during the annealing/extension step. Relative gene expression levels were quantified using the 2^−ΔΔCq method, with normalization to the internal reference gene 36B4 to minimize experimental variability (37).

Western blotting was performed according to established protocols (38). Briefly, cells were lysed using RIPA buffer (cat. no. SW104-02; Sevenbio) followed by sonication using a JY92-IIN ultrasonic cell disruptor (Ningbo Scientz Biotechnology Co., Ltd.) on ice at 20–25 kHz, 30% output power with cycles of 3 sec pulse and 3 sec interval, repeated for a total of 3 cycles to minimize protein denaturation. The resulting supernatants were collected, quantified using the BCA assay (Sevenbio), and denatured with 5× loading buffer (Beyotime Biotechnology). Protein samples containing 20 µg of total protein per lane were separated by SDS-PAGE using 7.5% polyacrylamide resolving gels. Following electrophoresis and transfer, PVDF membranes (MilliporeSigma) were blocked with 5% non-fat milk at room temperature for 2 h, then incubated with primary antibodies overnight at 4°C. Primary antibodies against YAP1 (1:10,000; cat no. 13584-1-AP) and β-actin (1:10,000; cat no. 66009-1-Ig) were purchased from Proteintech Group, Inc., with β-actin used as the loading control for normalization. After washing, the membranes were incubated for 2 h at room temperature with the appropriate HRP-conjugated secondary antibody (Proteintech Group, Inc.) diluted 1:100,000. Specifically, goat anti-mouse IgG (cat. no. SA00001-1) was used for mouse primary antibodies, while goat anti-rabbit IgG (cat. no. SA00001-2) was employed for rabbit primary antibodies. This was followed by protein band detection using an enhanced chemiluminescence detection kit from Beyotime Biotechnology and imaging with the ImageQuant 800 system (Cytiva).

Here, YAP1 ChIP-seq, H3K4me3, H3K27ac and ATAC-seq datasets were retrieved from the ChIP-Atlas database (chip-atlas.org/) to investigate the binding sites of YAP1 at the HRH1 promoter region. Based on the GENCODE database (GRCh38 assembly; http://www.gencodegenes.org/human/), the HRH1 RefSeq sequence is located on chromosome 3 at position 11,154,493-11,263,557 (https://genome.ucsc.edu/cgi-bin/hgSearch?search=HRH1&db=hg38). The downloaded BigWig files were visualized using the ‘Integrative Genomics Viewer’ (IGV; version 2.19.4; Broad Institute), with detailed dataset identifiers and download links provided in Table SIV. Since the ChIP-seq data presented in the present study were obtained exclusively from the ChIP-Atlas database, none of the cell lines were cultured or experimentally manipulated in the present study. Notably, the ‘PCa3’ in Table SIV refers to patient-derived prostate cancer organoids, which are well-characterized models derived from advanced prostate cancer metastases (39,40). Their genomic characterization (copy number alterations, mutations) and functional data are publicly accessible via the MSKCC cBioportal (http://www.cbioportal.org) and the NCBI PMC repository (pmc.ncbi.nlm.nih.gov/articles/PMC4237931/).

Pathway enrichment patterns were evaluated through GSEA implemented in the clusterProfiler R package (version 4.16.0; Bioconductor), analyzing transcriptomic data from TCGA-PAAD cohort (41). The Molecular Signatures Database was utilized to acquire predefined gene signatures, specifically the Cordenonsi_YAP_Conserved (M2871) and YAP1_UP (M2845) gene sets. Patients were stratified into high- and low-expression groups based on HRH1 levels and differentially expressed genes (DEGs) were identified using the ‘limma’ package (version 3.21, Bioconductor). The false discovery rate (FDR) was used for multiple-testing correction, with thresholds of FDR<0.05 and |log2FC|>1. Functional annotation of DEGs was conducted using Metascape (https://metascape.org/) to elucidate enriched biological pathways.

To predict chemosensitivity in TCGA-PAAD cohort, the present study employed the ‘oncoPredict’ R package (version 1.2; CRAN), a computational tool that integrates genomic features with drug response profiles (42). Training data were derived from the Cancer Cell Line Encyclopedia (CCLE; http://sites.broadinstitute.org/ccle/) and the Cancer Therapeutics Response Portal (CTRP; http://portals.broadinstitute.org/ctrp/), enabling the calculation of chemotherapy resistance scores for each patient.

The cellular composition of the TME was analyzed using the ‘deconvo_tme’ function within the ‘IOBR’ R package (43). This module integrates six well-established, publicly available deconvolution algorithms for the inference of cellular abundances from bulk transcriptomic data: CIBERSORT, MCP-counter, EPIC, xCell, quantTIseq and TIMER (44–49).

Furthermore, ssGSEA was performed to compute cellular enrichment scores for each patient using the ‘calculate_sig_score’ algorithm (50). Gene reference sets utilized for ssGSEA quantification were derived from authoritative literature (51–55). The TIDE computational method was employed to calculate the Immune Checkpoint Inhibitor (ICI) resistance score (55). Patients exhibiting high TIDE scores demonstrate susceptibility to immune evasion and exhibit lower response rates to immunotherapy.

To develop a robust and generalizable prognostic signature, the present study utilized the OmniLearn R package, which facilitates automated model construction across distinct machine learning algorithms (https://github.com/Feng-Rommel/OmniLearn). In line with previous reports (56), the present study implemented 101 specific algorithmic configurations (Table SV) to establish the prognostic model. These configurations were derived from 10 core survival modeling methods: Random Survival Forest (RSF), CoxBoost, stepwise Cox regression, Lasso, Ridge, Elastic Net (Enet), Survival Support Vector Machines (survival-SVM), Generalized Boosted Regression Modeling (GBM), Supervised Principal Components (SuperPC) and Partial Least Squares Cox regression (plsRcox). Models were established using single algorithms or paired combinations. In the combined strategy, genes were initially screened using the ‘RunML’ function (mode: ‘Variable’). Algorithms yielding <3 candidate features were excluded from further modeling, while those with ≥3 features proceeded to model construction using the ‘RunML’ function (mode: ‘Model’). Single-method models were constructed directly using the ‘RunML’ function (mode: ‘Model’).

Specific parameter settings for these ten core algorithms are described below. For stepwise Cox regression, all three directionality strategies (‘forward’, ‘backward’ and ‘both’) were tested. Penalized regression models (Lasso, Ridge, Enet) were fitted using the ‘glmnet’ package via the ‘cv.glmnet’ function, with the regularization parameter λ tuned through 10-fold cross-validation (CV). The elastic net mixing parameter α was systematically swept from 0 to 1 in increments of 0.1 (α=0: Ridge; α=1: Lasso; 0<α <1: Enet). Survival-SVM was implemented using the ‘survivalsvm’ package. GBM models were trained with the ‘gbm’ package under ten-fold CV. SuperPC, an extension of PCA adapted for survival outcomes, was executed via the ‘superpc’ package, with tuning performed using the ‘superpc.cv’ function (ten-fold CV). The plsRcox model was optimized using the ‘cv.plsRcox’ function from the ‘plsRcox’ package. Finally, RSF was implemented with the ‘randomForestSRC’ package using the ‘rfsrc’ function, fixing ntree to 500 and nodesize to 5 based on empirical validation to balance computational efficiency and predictive stability.

Following the construction of distinct prognostic models using TCGA-PAAD dataset (n=178) as the training set, risk scores for both the training cohort (TCGA-PAAD) and an external validation cohort (ICGC-PACA-AU; n=95) were calculated using the ‘CalPredictScore’ function. Subsequently, model performance was evaluated by computing the C-index via the ‘RunEval’ function. To ensure robustness, the optimal model was selected based on the highest mean C-index across both TCGA-PAAD and ICGC-PACA-AU cohorts. Finally, two microarray datasets (GSE28735 and GSE183795) were merged to constitute an additional external validation set (14,15). Time-dependent calibration was assessed using the Brier score at 1–4 years via the ‘pec’ package (Version 2025.06.24). This machine learning framework enables unbiased identification of robust prognostic biomarkers while mitigating overfitting through ensemble-based feature selection.

Uni- and multivariate Cox proportional hazards regression analyses, Kaplan-Meier (KM) curve generation and log-rank testing for survival differences were performed using the ‘survival’ (version 3.8–3; CRAN) and ‘survminer’ (version 0.5.1; CRAN) packages. Principal component analysis (PCA) was implemented via the ‘stats’ package and time-dependent receiver operating characteristic (ROC) evaluation was conducted using the ‘timeROC’ package (version 0.4; CRAN). For comparisons involving only two groups, unpaired Student's t-test or the Wilcoxon rank-sum test was used, as appropriate. When a single control group was compared against multiple treatment groups, statistical significance was assessed by one-way ANOVA followed by Dunnett's post hoc test. The correlation between two variables was investigated using a non-parametric Spearman test. All statistical tests are two-tailed. Data from experiments are presented as mean ± standard deviation (SD). All statistical analyses were performed in R (version 4.5.1), and P<0.05 was considered to indicate a statistically significant difference.

To identify potential therapeutic targets for PDAC, the present study intersected three gene sets: 381 GPCR genes, 1,285 genes significantly associated with prognosis (log-rank P<0.01) in TCGA-PAAD cohort and 2,456 genes markedly upregulated (log2FC >2) in PDAC from the same cohort. This analysis yielded two candidate genes, GPRC5A and HRH1 (Fig. 1A). Due to the lack of specific pharmacological agents targeting GPRC5A, HRH1 was chosen as the main focus of this investigation because it has targeted inhibitors (57).

KM analysis in TCGA-PAAD cohort revealed that high expression of HRH1 was significantly associated with unfavorable overall survival (Fig. 1B). Both univariate and multivariate Cox regression analyses further identified HRH1 as an independent prognostic factor in pancreatic cancer (Fig. S1). Analysis using the GEPIA2 online platform confirmed significant overexpression of HRH1 in PDAC samples (Fig. 1C). Moreover, bulk RNA sequencing of 266 matched patient samples across seven independent cohorts demonstrated consistently higher HRH1 expression in tumor tissues relative to adjacent non-tumor counterparts (Fig. 1D). Analysis of the GSE26088 dataset revealed that HRH1 expression was elevated across 19 pancreatic cancer cell lines relative to the normal pancreatic ductal epithelial cell line HPDE (Fig. 1E). qPCR validated increased HRH1 expression in seven PDAC cell lines compared with the non-cancerous HPNE cell line (Fig. 1F). At single-cell resolution, malignant cells exhibited higher average UMI counts for HRH1 than normal ductal cells (Fig. 1H).

TMB analysis suggested a positive association between high HRH1 expression and increased mutation load (Fig. 1G). Evaluation using the ‘TIDE’ algorithm revealed a positive association between HRH1 expression and TIDE score, implying that high HRH1 expression may be associated with diminished ICI efficacy (Fig. 1I). Additionally, chemoresistance prediction analysis indicated that high HRH1 expression (stratified by median expression) was associated with resistance to multiple chemotherapeutic agents, including gemcitabine, fluorouracil, SN-38, oxaliplatin, doxorubicin, platinum-based drugs, mitomycin and vincristine (Figs. 1J-M, S2).

To further investigate the potential therapeutic relevance of HRH1 in PDAC, the present study performed MR analysis. The IVW method indicated that fexofenadine treatment is associated with a reduced risk of pancreatic cancer (OR=1.71×10⁻²⁰; Figs. 2A, 3A; Table SVI). Leave-one-out sensitivity analysis confirmed the robustness of these findings (Fig. 2B). The analysis showed no significant heterogeneity (MR-Egger Q=9.224) and no evidence of horizontal pleiotropy (Fig. S3B; Table SVII).

Following HRH1 knockdown using siRNA, the proliferation of PANC-1 and SW1990 cells was inhibited compared with control cells (Fig. 2C and D). In PANC1 cells, silencing HRH1 enhanced sensitivity to gemcitabine (Fig. 2E). The IC₅₀ for the siControl group was 163.9 µM, whereas it was markedly reduced to 8.226 µM and 3.645 µM in the siHRH1#1 and siHRH1#2 groups, respectively (Fig. 2E). A similar potentiation of gemcitabine cytotoxicity was observed in SW1990 cells upon HRH1 knockdown. The IC₅₀ value for the siControl group was 3.827 µM, which decreased to 0.2241 µM and 0.09096 µM in the siHRH1#1 and siHRH1#2 groups, respectively (Fig. 2F). The RT-qPCR validation data for HRH1 knockdown efficiency are provided in Fig. 3E and F. Collectively, these results underscore HRH1 as a promising therapeutic target in pancreatic cancer.

To investigate the signaling pathways activated by HRH1 in PDAC, the present study divided samples into high and low expression groups based on the median expression of HRH1 and identified DEGs between the two groups. Metascape enrichment analysis revealed that upregulated DEGs in the high HRH1 expression group were associated with YAP1 signaling (Fig. 3A). GSEA analysis indicated a positive correlation between HRH1 and the ‘Cordenosi YAP conserved signature’ as well as the ‘YAP1 up signature’ (Fig. 3B and D). In TCGA-PAAD cohort, HRH1 expression showed a positive correlation with the expression of multiple target genes of the YAP1 signaling pathway (Fig. 3C). Knockdown of HRH1 in human endothelial cells downregulated multiple YAP1 target genes, including CCN1 and CCN2 (Fig. S4A). These results indicated that HRH1 promotes activation of the YAP1 signaling pathway in endothelial cells. In PANC1 and SW1990 cells, knockdown of HRH1 resulted in decreased expression of CTGF, CYR61 and ANKRD1, while the transcriptional levels of YAP1 remained largely unchanged (Fig. 3E and F), indicating that HRH1 may influence YAP1 pathway activity by affecting YAP1 protein stability rather than its transcription in pancreatic cancer cells. Furthermore, in breast cancer (MDA-MB-231; GSE59232), neuroblastoma (NLF; GSE130401), liver cancer (HepG2; GSE137915), and renal clear cell carcinoma (RCC4; GSE146354), silencing YAP1 led to a significant decrease in HRH1 expression at the transcriptome level (Fig. S4B-E). Similarly, through qPCR experiments, it was confirmed that knockdown of YAP1 in PANC1 and SW1990 cells led to decreased expression of HRH1 mRNA (Fig. 3E and F).

YAP1 is a transcriptional co-activator that binds to TEA domain transcription factors (TEAD) proteins and anchors to the promoter regions of target genes to regulate downstream expression (58). Further analysis using ChIP-Atlas data indicated enrichment of H3K4me3 (an active promoter marker) and H3K27ac (an enhancer/promoter marker) in the genomic region ~1,000 bp upstream of the HRH1 transcription start site (59,60). This region also exhibited high signal in ATAC-seq (open chromatin), suggesting it is likely the promoter region of HRH1 (Fig. 4A). YAP1 ChIP-seq data revealed peaks in the HRH1 promoter region, indicating direct binding of YAP1 to the HRH1 promoter (Fig. 4A). In PA-TU-8902 (pancreatic cancer cells) and PCa3 (prostate cancer cells), knockdown of YAP1 and TAZ led to reduced chromatin accessibility in the HRH1 promoter region (Fig. 4B). Notably, although not all cell lines in Fig. 4A and B are pancreatic cancer-derived, YAP1 binding to the HRH1 promoter was consistently detected across these diverse cell types, indicating that this regulatory mechanism is broadly active across multiple human cancer types and not restricted to pancreatic cancer cells.

GPCR activation transmits signals through the Rho/ROCK/F-actin pathway, thereby regulating YAP1 activity (61). Western blot analysis showed that knockdown of HRH1 in PANC1 and SW1990 cells decreased YAP1 protein expression levels (Fig. 4C and D; Table SVIII). Moreover, inhibition of HRH1 with fexofenadine decreased YAP1 expression compared with the control group, while activation of histamine receptors with histamine increased YAP1 expression. Co-treatment with fexofenadine and histamine also led to decreased YAP1 expression, though to a lesser extent than with fexofenadine alone (Fig. 4E and F; Table SIX). Taken together, these findings indicate that HRH1 stabilizes and activates YAP1 signaling, whereas YAP1 knockdown downregulates HRH1 transcription, potentially through direct binding to the HRH1 promoter region.

Given the key role of the HRH1/YAP1 signaling axis in pancreatic cancer, the present study employed 101 distinct machine learning algorithms to construct prognostic models based on HRH1, YAP1 and YAP1 downstream genes associated with HRH1 expression (listed in Fig. 3D). A prognostic model was first constructed based on TCGA-PAAD training cohort, and individual risk scores were then generated for patients in both TCGA-PAAD and ICGC-PACA-AU cohorts, with the latter serving as an independent validation set. Batch effects were removed across different datasets (Fig. S5). The C-index was employed as the main evaluation metric to determine the model's predictive performance, reflecting the likelihood of agreement between model predictions and actual observed outcomes. Among the 101 algorithms evaluated, 25 models were successfully constructed. In TCGA cohort, the mean C-index was 0.659 (SD=0.069), while in the ICGC cohort, the mean C-index was 0.683 (SD=0.08; Fig. 5A). Among these, the ‘plsRcox’ approach, which combines partial least squares regression with Cox regression, achieved optimal predictive performance, exhibiting a mean C-index of 0.713. The HRH1/YAP1 signaling-derived prognostic signature comprises eight genes. Risk scores were computed based on a weighted linear combination of expression levels from eight signature genes, defined as: risk score=(0.032 × HRH1) + (0.015 × YAP1) + (0.227 × ECT2) + (0.215 × ITGB5) + (0.207 × SHCBP1) + (0.108 × MAML2) + (0.032 × YWHAZ) + (0.147 × ITGA2). Analysis revealed that the risk score based on HRH1/YAP1 signaling emerged as an independent determinant of overall survival in pancreatic cancer patients, as shown by univariate and multivariate Cox regression (Fig. 5B and C).

To further validate the robustness of the model, two additional GEO datasets (GSE28735 and GSE183795), which were not used during model training, were employed as external validation sets. Higher risk scores were associated with increased mortality (Fig. 6A). Using the median risk score of 6.171 from TCGA-PAAD training cohort as the cut-off, both the training and validation sets were stratified into high- and low-risk groups. KM analysis revealed significantly worse survival in the high-risk group, with log-rank P<0.001 in TCGA-PAAD and GEO cohorts, and P<0.05 in the ICGC-PACA-AU cohort (Fig. 6B). Time-dependent ROC analysis demonstrated that the AUC values for predicting 1-, 2-, 3- and 4-year survival were 0.778/0.734/0.699/0.699 in TCGA-PAAD, 0.786/0.835/0.570/NA in ICGC-PACA-AU and 0.593/0.688/0.813/0.813 in GEO, respectively (Fig. 6C). The time-dependent Brier scores for 1-, 2-, 3- and 4-year overall survival prediction were 0.143, 0.185, 0.190, and 0.182 in TCGA-PAAD cohort; 0.115, 0.185, 0.202, and 0.182 in the ICGC-PACA-AU cohort; and 0.242, 0.187, 0.128, and 0.100 in the GEO cohort, respectively. Lower Brier scores indicate improved calibration and predictive accuracy (Fig. S6). Heatmap analysis demonstrated a positive correlation between all model genes and the risk score (Fig. 6D). PCA analysis showed clear separation between high- and low-risk groups, with non-overlapping distributions in the scatter plot (Fig. 6E).

Previous studies have reported that targeting HRH1 can enhance the sensitivity to ICIs by upregulating MHC-I expression in pancreatic cancer, while targeting YAP1 increases chemosensitivity to gemcitabine (62,63). Inspired by these findings, the present study further investigated the association between the HRH1/YAP1 signaling-related risk score and sensitivity to both chemotherapeutic agents and ICIs. As illustrated in Fig. 7A, high-risk patients demonstrated markedly elevated IC₅₀ levels for standard chemotherapeutic agents used in PDAC treatment. These included gemcitabine, fluorouracil, SN-38 (the active derivative of irinotecan), and oxaliplatin, along with other commonly administered drugs such as platinum-based compounds, doxorubicin, mitomycin and vincristine. Further analysis revealed a higher mutation frequency of the pancreatic cancer driver gene KRAS in the high-risk group (85 vs. 44%), with similar trends observed for TP53 (78 vs. 44%), SMAD4 (26 vs. 16%) and CDKN2A (26 vs. 11%). The heightened chemoresistance and worse prognosis evident in high-risk patients may be partially attributable to mutations within these genes (Fig. 7B) (5,64).

Although TMB is generally associated with improved response to immunotherapy, with a typical cut-off for TMB-high defined as 10 muts/Mb (65–67). The present analysis revealed a positive correlation between the risk score and TMB. However, the majority of samples demonstrated a TMB below 4 muts/Mb. This suggests that these patients are unlikely to respond favorably to ICIs (Fig. 7D). The present analysis of the TME further revealed that a high risk score was associated with an enriched presence of immunosuppressive cells, including myeloid-derived suppressor cells (MDSCs), cancer-associated fibroblasts (CAFs), regulatory T cells (Tregs) and T helper 2 cells (Th2 cells) (Figs. 7C and S7). In contrast, this high-risk phenotype was characterized by a significant decrease in CD8⁺ T cells. The TIDE algorithm, designed to forecast responses to ICI treatment, yielded elevated scores in the high-risk group (Fig. 7E). Elevated TIDE values are associated with a greater likelihood of unfavorable therapeutic outcomes. Therefore, although pancreatic cancer is generally considered an ‘immune-cold’ tumor with limited response to ICIs, the low-risk group benefiting from immunotherapy.