Samples were collected in several sites along the country: Clínica Indisa, Biobanco de Tejidos y Fluidos de la Universidad de Chile, Hospital Padre Hurtado, Hospital Regional de Concepción, Hospital Regional de Talca, Hospital de Puerto Montt, Hospital San Juan de Dios, Instituto Nacional del Cáncer, Hospital Regional de Coquimbo, and Hospital Regional de Arica. Hematoxylin and eosin-stained tissue sections, derived from Formalin-Fixed Paraffin-Embedded (FFPE), underwent a thorough review by two pathologists. The tumor region was delineated. Only specimens with a percentage greater than 10% of tumor content were eligible for inclusion in this study. Areas characterized by significant necrosis and intra-tumoral fibrosis were excluded from the analysis. Unstained tissue sections (5 µm-thick) were prepared, deparaffinized, and the tumor tissue specimen was carefully collected for subsequent DNA extraction.

DNA was extracted from paraffin-embedded samples using the RecoverAll™ Total Nucleic Acid Isolation Kit (Thermo Fisher Scientific). A 20 ng DNA input was used for library preparation with the Oncomine™ Comprehensive Assay V1 (OCAv1, Thermo Fisher Scientific), following the manufacturer’s protocol. Library preparation was performed using the Ion Chef™ System, with a final library concentration of 100 pM. Purified DNA was quantified using the Qubit™ dsDNA HS Assay and the Quant-iT™ PicoGreen^® dsDNA Reagent Kit (Invitrogen, Thermo Fisher Scientific, Carlsbad, CA, USA). DNA purity was evaluated by measuring the absorbance ratio at 260/280 nm. The DNA integrity and fragmentation status were assessed using the High Sensitivity Genomic DNA Analysis Kit (DNF-488) on a Fragment Analyzer system (Agilent Technologies, formerly Advanced Analytical). Sequencing was conducted on an Ion PGM™ System in single-read mode using the Ion 550 kit-CHEF (Thermo Fisher Scientific). The list of genes analyzed and the sequencing quality metrics for all analyzed samples, including total reads, coverage, and uniformity, are summarized in Supplementary Tables S1 and 2 . Of the 118 FFPE tumor samples collected, 66 (55.9%) yielded sufficient quantity and quality of DNA to proceed to sequencing. The remaining 52 samples (44.1%) failed pre-sequencing QC due to low DNA yield, high degradation. Many of the samples that failed due to low DNA yield or high degradation were obtained from various pathology centers. In these cases, inconsistencies in pre-analytical processes, such as prolonged fixation times and the use of non-buffered formalin, may have significantly contributed to DNA degradation. Ten samples were excluded during post-sequencing quality control for failing to meet predefined variant calling thresholds (≥3 million total reads and ≥90% on-target rate). The final genomic analysis was therefore performed on the remaining 56 patient samples ( Supplementary Figure S1 ).

The data preprocessing and processing were carried out using the OCAv1 v5.18 DNA workflow, using default parameters of the GRCh37/hg19 reference genome. For alignment and variant calling, stringent parameters were defined. Single Nucleotide Variants (SNV) and Indels required a minimum allele frequency of 0.05 and 0.07. The minimum coverage for a variant to be considered was set at 10x for SNV and Indels. Additionally, the minimum coverage for the variant location was set at 200x. Variant annotation was performed using ANNOVAR (9), including RefGene, GnomAD v2.1.1, ESP6500, 1000 Genomes phase 3, CADD v1.3, dbSNP v150, COSMIC v94, CLINVAR 2021, ICGC28, ABraOM, and Snp138NonFlaged. To enhance the filtering of germline variants in tumor samples, large and local population germline variant databases were interrogated:CSVS (10), GnomAD (overall and population specific); BIPMed (Brazilian Initiative on Precision Medicine); and a Chilean database (3, 11) (variants imputed from genotyping two arrays with 1,313 and 2,249 samples from Chilean individuals, and whole exome sequencing (WES) data from 87 individuals (NCBI dbSNP database under accession code 1062069) (12).

Mutation data of GBC samples were extracted from the Memorial Sloan Kattering Clinical Sequencing cohorts (6) through cBioportal (https://www.cbioportal.org, accessed on march 1st 2024). In addition, GBC mutation data from Japanese and Singaporean samples were extracted from the International Cancer Genome Consortium (ICGC) website (downloaded on march 1st 2024). Only somatic protein-affecting variants found in primary tumors were analyzed in all datasets.

The mutational landscape of GBC samples was characterized using the R/maftools package (13). To identify actionable and driver variants within the cancer samples, the Cancer Genome Interpreter platform (https://www.cancergenomeinterpreter.org/home) was employed. Variants were labeled according to the level of evidence (A, B, C, and D) obtained from VICC integrated knowledge base, following the AMP/ASCO/CAP guidelines. OncodriveMUT and BoostDM algorithms were used for prediction of driver variants (14). Only predicted drivers that were also annotated as “oncogenic or likely-oncogenic” in OncoKB (oncokb.org) were considered.

Genetic ancestry was evaluated using genomic DNA extracted from peripheral blood samples of 37 individuals. The ancestry analysis was performed using the ADMIXTURE software (15) and Genetic PCA with the Eigenstrat function (16). This approach allows for supervised estimation of individual ancestry components, including African, European, Native American Mapuche and Aymara, using a large number of Single Nucleotide Polymorphism (SNP) markers previously defined as Ancestry informative marker (AIM) (15). The panel for the preselection of ancestry-informative markers and the estimation of individual ancestry proportions in genetically admixed Chileans, which included 63 Aymara individuals (17) and 28 Mapuche individuals (17, 18). Additionally, it comprised 206 Europeans (99 Utah residents with Northern and Western European ancestry [CEU] and 107 Iberians from Spain [IBS]) and 108 African Yoruba from Ibadan, Nigeria (YRI) from the 1000 Genomes Project (17).

A Chi-square test was performed to compare mutation frequencies between the Chilean GBC cohort and those from Japan (GBC-JPN), Singapore (GBC-SGN), and the MSK2022 cohort. p-values were adjusted using the Benjamini–Hochberg method to control the false discovery rate (FDR), given the large number of genes tested simultaneously. Pairwise Chi-square tests with Yates’ continuity correction were performed as post hoc comparisons in cases where the adjusted p-value (FDR) was below 0.05. Multivariate logistic regression models were adjusted to evaluate the association between ancestry and the presence of variants in target genes. Sex, age at the time of collection, history of cholelithiasis, and GBC histology were included as covariates using the following model:

In this model, each target gene status (mutated/non-mutated) was modeled as a function of each ancestry, grouped as high/low, along with the effect of other covariates mentioned. A significance level of 0.05 was defined as statistically significant. Patients were categorized into two groups using the median of each estimated genetic ancestry (European, Mapuche, Aymara, and African). All statistical analyses were conducted using R 4.2.1 software.

We collected a total of 118 FFPE biopsies from GBC. To account for the high heterogeneity of the Chilean population, samples were gathered from several medical centers throughout the country. Of these, 94 (76.6%) were from female patients and 24 (23.4%) from male patients, respectively, with an average age of 62 years at the time of collection. Seventy-eight patients (66.1%) had a history of cholelithiasis in their clinical history ( Table 1 ).

In terms of ethnicity, most patients self-declared Chilean (Admixed) (79.8%), followed by Mapuche (8.5%), Aymara (3.4%) and European (3.4%). Genetic ancestry distribution at individual level is provided in Supplementary Figure S3 . The specific genetic ancestry As for the histological diagnosis of the samples, a varied distribution in this type of cancer was revealed, with adenocarcinoma being the most frequent type (84.9%), while squamous cell carcinomas and high-grade dysplasia were less common. Most patients were diagnosed at G1 and G2 stages at the time of surgery ( Table 1 ).

Just 56 of the 118 samples met all quality control criteria for DNA, library, and sequencing data ( Supplementary Figure S1 ). This finding underscores the significant influence of tissue processing and other pre-analytical factors on the successful execution and quality of molecular analyses. We identified 535 somatic mutations across 43 tumor samples. The majority of these somatic variants, 469 (87.66%), were missense mutations, followed by 34 (6.36%) nonsense mutations, 12 (2.24%) splice site mutations, and 5 (0.93%) frameshift deletion mutations. ( Figure 1A ), with the majority falling under the Single Nucleotide Variants (SNV) category ( Figure 1B ). The C > T nucleotide variation was the most frequent substitution among the detected somatic variants. Additionally, transitions (Ti) are more prevalent than transversions (Tv) among the SNVs identified in all genes with mutations ( Figures 1C, D ).

Somatic mutations were identified in 41 out of the 161 genes analyzed. The most frequently mutated genes included TP53 (30%), TSC2 (29%), NOTCH1 (27%), BRCA2 (16%), KRAS (16%), and PTCH1 (16%) ( Figure 2 ). The customized oncoplot integrates both gene-level and pathway-level alterations, grouping genes according to their corresponding oncogenic signaling pathways. In this analysis, 2/6 (33.3%) genes in the TP53 pathway and 15/85 (17.7%) genes in the RTK–RAS pathway were altered. Notably, the RTK–RAS pathway was affected in 50% of GBC samples, followed by the TP53 pathway (38%), NOTCH (37%), and PI3K (32%) ( Figure 2 ). Among the 10 most mutated genes, several genes show a significant co-occurrence, such as MLH1/TSC2, NOTCH1/MLH1, and NOTCH1/TSC2, among others. Mutations in TP53 and KRAS seem to be mutually exclusive, although this trend did not reach statistical significance ( Figure 3 ).

Given the limited therapeutic options for patients with GBC, identifying targets for therapies that have already demonstrated clinical benefit in other tumor types may enable the inclusion of GBC patients in ongoing or future trials. In our cohort, all 44 patients with somatic variants harbored one or more actionable alteration with supporting Level A evidence. The most frequently altered genes were KRAS wildtype (n=35, predictive of response to anti-EGFR therapies such as panitumumab and cetuximab in colorectal cancer), TSC2 (n=16, associated with everolimus in renal angiomyolipoma and giant cell astrocytoma), and BRCA2 (n=9, associated with PARP inhibitors such as rucaparib in ovarian, prostate, and pancreas cancers). Additional Level A alterations were identified in ATM, BRCA1, EGFR, ERBB2, FGFR2, FGFR3, PTCH1, RET, and TSC1, all of which predict responsiveness to specific targeted therapies. Importantly, we detected KRAS mutations (G12D/V/C, G13D, Q61H) associated with resistance to EGFR inhibitors such as cetuximab and panitumumab in colorectal and lung cancer, and responsiveness to KRAS inhibitors sotorasib and adagrasib in colorectal, NSCLC, pancreatic, ampullary, and hepatobiliary cancer ( Table 2 ).

Although genomic data on GBC are available in public repositories, this data is still scarce and has a limited representation of Latin American patients. Given the marked geographical disparities in GBC incidence and mortality, it is plausible to hypothesize that underlying environmental and genomic background may favor the activation of distinct molecular oncogenic mechanisms in high incidence regions. To explore this possibility, first we compared the mutation frequencies in the 41 genes identified in the Chilean cohort with those reported in previously published datasets MSK-IMPACT 2022 cohort, as well as two GBC cohorts from Japan and Singapore (obtained from the International Cancer Genome Consortium (ICGC) (https://dcc.icgc.org/). Despite the limited sample size, we found significant differences in the mutation frequencies of all genes but ERBB2 across cohorts ( Table 3 ).

Next, we identified variants predicted and annotated as oncogenic and likely-oncogenic in 23 (41.07%) out of the 56 Chilean samples analyzed. These mutations are distributed in 12 genes. TP53 and KRAS exhibit the highest frequencies of oncogenic/likely-oncogenic variants in 18% and 9% of the studied samples, respectively, followed by ERBB2 (5%) ( Supplementary Figure S2 ). TP53 variants are predominantly missense mutations situated within the p53 DNA-binding domain ( Figure 4A ). KRAS recurrent oncogenic G12D/C/V, G13D, and Q61H mutations were found ( Figure 4B ). In ERBB2 gene, activating S310F, R678Q, and L755S were identified at Furin-like, TMD, and Kinase domains ( Figure 4C ).

To further explore whether the overall differences detected in gene mutation frequencies between GBC from mostly White (MSK.,2022), Asian (GBC-JPN, GBC-SGNAs), and Chilean patients involved key driver mutations, oncogenic variants in the 12 driver genes were identified and compared in all cohorts. Most oncogenic mutations were also found in one or more studies, especially hotspot mutations in CDKN2A, ERBB2, SMAD4, FBXW7, KRAS, and TP53. Nevertheless, differences are also evident, in tumor suppressor genes ATM and TP53. Notably, in Chilean samples, likely oncogenic EGFR R776C, BRCA2 E3002K, and MLH1 W666* and Q562*, were found. No oncogenic mutations were identified in these genes in the referenced studies ( Table 4 ).

Genomic studies have provided evidence supporting that Chilean Amerindian genetic ancestry is a risk factor for GBC (3). In fact, a correlation has been established indicating that Chileans with a higher proportion of South Amerindian genetic ancestry (Mapuche, MAP) have a higher risk for GBC (3, 8). With the aim to assess whether this genetic-ancestry-related risk for GBC is also associated with the risk for a specific molecular oncogenic mechanism of disease, we used a panel of Ancestry informative marker (AIM) to estimate the proportion of ancestral genetic backgrounds for 37 out of 43 patients with somatic variants. As expected, European and Mapuche genetic ancestries are predominant in the analyzed individuals (0.493, sd: 0.145; and 0.34, sd: 0,129, respectively), followed by Chilean main North Amerindian Aymara (0.137; sd: 0,221), and African (0.019; sd: 0.16) ancestries ( Supplementary Figure S3 ).

A general overview of the somatic mutation´s data suggest that higher proportion of Mapuche ancestry would associate with a lower number of somatic mutations ( Supplementary Figures S4 , S5 ). Although these results are preliminary, this led us to investigate the potential association between Mapuche ancestry and specific gene mutation frequency in TP53 (the most frequently mutated gene in our restricted cohort). Higher proportion of Mapuche ancestry may be positively associated with mutations in TP53; however, the result is not significant (OR = 6.73, 95% CI: 1.14, 59.6, p = 0.053). In addition, poorly differentiate (G3) tumors showed positive associations regarding the mutational status of TP53 (p= 0.006) in a multivariate analysis ( Table 5 ).

Beyond our findings, studies in other high-incidence regions, such as by Mishra et al. in India, similarly show a high frequency of TP53 mutations (90.90%), alongside other common alterations like SMAD4, NOTCH1, and ERBB2 (22). This shared high frequency of TP53 mutations between Chilean and Indian GBC patients may stem from common risk factors like gender, bowel habits, and cholelithiasis (23). Chronic inflammation induced by cholelithiasis is a plausible mechanism, potentially triggering TP53 mutations that contribute to GBC development (24). While a common variant in the ABCB4 gene (rs4148808) linked to gallstone disease might play a role in both populations, further studies are needed for verification (25). It is crucial to note that differences in variant allele frequency (VAF) thresholds (e.g., >1% in Mishra et al. vs. >5% in our study) could influence observed mutation frequencies and interpretations of somatic variants.

We also observed a high number of non-synonymous somatic mutations in two patients, both carrying variants in DNA repair genes such as MLH1 and MSH2. Alterations in these genes have been associated with increased mutational load and may suggest potential sensitivity to immune checkpoint inhibitor (ICI) therapies (26). Furthermore, the predominance of C > T substitutions in our mutation pattern indicates a possible APOBEC-mediated mutational signature (27), a finding also reported in a whole exome sequencing study including Chilean GBC patients (21). This opens avenues for signature-based therapies.

Our cohort showed a female-to-male ratio of 4:1, consistent with the well-established female predominance in GBC incidence reported in Chile and other high-incidence regions (1–3). This disparity is thought to be largely driven by the higher prevalence of gallstones and chronic cholecystitis among women, potentially influenced by hormonal and metabolic factors (1–3). While sex-related differences in mutational patterns were not statistically evaluated in our series due to sample size constraints, previous studies have suggested that certain genomic alterations in GBC may vary by sex, including recent evidence describing distinct clinical and molecular profiles between men and women (28). However, these differences are often confounded by underlying risk factor distributions rather than direct sex-linked pathogenic mechanisms. Nonetheless, the integration of sex as a biological variable in future genomic studies of GBC could help clarify its role in tumor development and progression.

Given the limited approved targeted therapies for GBC and the reliance on systemic chemotherapy, we explored the actionable potential of genetic alterations in our Chilean cohort. Using the CGI platform, we identified clinically actionable alterations with drug recommendations in genes such as ATM, BRCA1/2, EGFR, ERBB2, FGFR3, KRAS, PTCH1, RET, and TSC1/2. While many of these drugs are established in other solid tumors, their specific utility in GBC is emerging. For instance, panitumumab combined with gemcitabine and oxaliplatin showed promising results in KRAS wild-type biliary tract cancer (29), and the SGNTUC-019 study demonstrated the efficacy of Tucatinib and Trastuzumab in HER2-positive metastatic biliary tract cancer (30). The recent FDA approval of sotorasib, a KRAS G12C inhibitor, further highlights the potential for targeted therapies (31).

Notably, ATM, BRCA1, and BRCA2 mutations are highly significant due to their involvement in the DNA damage response (DDR) pathway. PARP inhibitors (e.g., olaparib, niraparib) have shown efficacy in BRCA-mutated cancers (32), with recent evidence suggesting benefit in GBC patients with BRCA2 alterations (33). Although ATM alterations are less studied in GBC, their reported presence in biliary tract tumors suggests they could serve as predictive biomarkers (34, 35). These findings collectively underscore the potential for implementing personalized therapies for Chilean GBC patients by repurposing existing drugs based on specific genomic alterations.

The recurrent involvement of RTK-RAS, PI3K, and TP53 pathways in our cohort underscores their central role in the biology of BTC. Alterations in the RTK-RAS pathway, including KRAS mutations, are known to drive tumor proliferation and survival, and may influence sensitivity to targeted agents or MEK inhibitors (36). Similarly, PI3K pathway activation, through mutations in PIK3CA or upstream receptor tyrosine kinases, has been implicated in resistance to standard chemotherapy, and PI3K/AKT/mTOR inhibitors are currently being evaluated in BTC clinical trials (37). TP53 inactivation, the most frequent event in our cohort, is not only associated with loss of cell cycle control and genomic instability but has also been linked to reduced sensitivity to gemcitabine-based chemotherapy in BTC. This highlights p53 as a potential therapeutic target, with different strategies depending on mutational status: in TP53-mutated BTC, direct reactivation of mutant p53 or synthetic lethality approaches targeting cell cycle checkpoint proteins (e.g., Chk1, ATR, Wee1) may be promising; in wild-type p53 BTC, inhibition of its negative regulators (e.g., MDM2, WIP1) could enhance p53 activity (38). Such pathway-level characterization in GBC not only improves our understanding of tumor biology but also informs the rational design of targeted and combination therapies, particularly in the context of overcoming chemoresistance.

Using computational algorithms, we identified driver variants in TP53, KRAS, and ERBB2, genes frequently mutated in other biliary tract tumors (39). TP53 variants, often influenced by external mutagens, include R280 (associated with aristolochic acid exposure) (40, 41) and R248 (linked to tobacco smoke in lung cancer) (42). The recurrence of R248W and R282W variants may stem from increased CpG site methylation during inflammation, a process contributing to GBC tumorigenesis (43, 44).

Our comparison with international cohorts revealed that ERBB2 was the only gene without significant mutation frequency differences across our Chilean cohort and those from the USA, Singapore, and Japan (45, 46). This consistent frequency, also reported by Nepal et al. (21) and Bitter et al. (47), and, importantly, the frequency of ERBB2 alterations was also found to be comparable between Chilean and U.S. cohorts (48). It has been reported that GBC patients harboring ERBB2 alterations exhibit improved overall survival. Thus, ERBB2 may represent a broadly applicable prognostic and therapeutic biomarker across diverse populations (48).

Beyond common drivers, our study uncovered distinct patterns. Two patients exhibited hypermutated tumors, potentially linked to variants in DNA repair genes (MLH1, MSH2) and suggesting responsiveness to immune checkpoint inhibitor therapies (26). The predominance of C > T substitutions in our cohort also points to an APOBEC-mediated mutational signature (27), consistent with previous findings in Chilean GBC patients (21), opening avenues for signature-based therapeutic strategies.

Furthermore, we explored the association between Mapuche ancestry and TP53 mutation likelihood. While not statistically significant, our exploratory analysis suggested a possible association between higher Mapuche ancestry and TP53 mutation status (OR = 6.73, p = 0.053). While this trend did not reach statistical significance, it may reflect underlying biological differences warranting further investigation. This observation should be considered hypothesis-generating, and validation in larger, independent GBC cohorts will be necessary to confirm or refute this potential ancestry-related effect. This aligns with prior research linking Mapuche ancestry to GBC risk variants (e.g., in ABCG8, TRAF3, ABCB1, ABCB4 genes) and increased susceptibility to gallstones, a known inducer of TP53 mutations via chronic inflammation (24, 25, 49, 50). Such observations underscore the influence of genetic ancestry on tumor profiles, as seen in pan-cancer analyses (51).

These results suggest that there may be tumor profiles associated with genetic ancestry. Although the results are preliminary, they indicate that patients with greater Mapuche ancestry may have a lower mutation rate, which could reflect biological differences in the tumorigenesis of this type of malignant tumor, supporting the need to incorporate genomic studies in GBC.

Recently, Zhu et al., reported overall similar mutation profiles in a Chilean vs US GBC cohorts. However, the immune profiles were highly distinct, characterized by higher densities of T cells and PD-1, but lower macrophages in tumors from Chilean patients (52). Although these clinically relevant differences may reflect lifestyle and/or environmental factor dissimilarities between the 2 countries, their findings reinforce the idea that the genetic origin of patients could modulate cancer biology and response to treatment, thus supporting the development of precision medicine strategies that incorporate ancestry, especially in underrepresented populations such as Chileans.

Given the scarcity of approved GBC therapies, identification of actionable alterations (e.g., ATM, BRCA1/2, EGFR, ERBB2, FGFR3, KRAS, PTCH1, RET, TSC1/2) in Chilean patients provide a rationale for personalized treatment, particularly as specific therapies for KRAS and HER2-positive BTC are emerging (29–31). Mutations in ATM and BRCA1/2, crucial for DNA damage response, highlight the potential for PARP inhibitors, with promising results already seen in a GBC patient with a BRCA2 mutation (32, 33).

Despite these insights, our study has limitations. For instance, the NGS assay used (Oncomine Comprehensive Assay v1 (OCAv1), although it covers a panel of 161 clinically relevant cancer genes, it does not detect alterations in, for example, ARID1A, which is frequently mutated in GBC, including in Chileans (21, 52, 53).

In addition, although this assay is designed to detect a certain number of structural variants, only a few samples passed QC for RNA sequencing. Same with CNV. This could lead to underestimation of clinically relevant events such as ERBB2 amplification (49). The future application of broader sequencing approaches, such as whole-exome sequencing (WES) or whole-genome sequencing (WGS), in Chilean GBC cohorts would allow for a more comprehensive characterization of the mutational spectrum, the discovery of novel prognostic or therapeutic biomarkers, and the strengthening of analyses by integrating the ancestral component.

Nevertheless, challenges with FFPE sample quality have to be urgently addressed to further advance in a comprehensive and reliable genomic characterization of GBC in the country. Optimizing biopsy handling and preanalytical conditions are crucial to minimize sample attrition and enhance genomic profiling (54–56). Future work should also consider liquid biopsies (cfDNA) as a complementary strategy, given their concordance with tissue mutations and potential for non-invasive detection of actionable variants (57).

In conclusion, we identified actionable and driver genetic variants that may inform tumor response in Chilean GBC patients, laying the groundwork for personalized GBC therapies. Future research should leverage comprehensive genomic methodologies, such as whole-exome sequencing (WES), whole-genome sequencing (WGS), and liquid biopsies, to overcome current limitations and further elucidate GBC’s complex genetic and environmental interactions.