High-throughput drug screening

A total of 18 384-well micro-well plates were used in this study, two for each organoid line, on which 11 drug combinations were dispensed in a 4 × 6 drug synergy matrix (AF x DrugX).

The growth rate of different organoid lines was measured over 120 h using image-based quantification. It was clear that baseline growth rate varied significantly across the different organoid lines (Fig. 1A), thus highlighting its importance when evaluating the effects of drugs and drug combinations on organoids. Therefore, we employed a growth rate-based metric to evaluate therapeutic effect, the normalized-organoid growth rate (NOGR), which considers basal organoid growth rate as well as viability (Fig. S1, Additional file 1). At day five, 5472 unique NOGR datapoints were obtained for all 9 lines of which 147 (2.69%) were identified as outliers in relation to the fitted drug-response curve and excluded from further analysis. A high correlation was observed between the non-fitted and fitted NOGR values (r = 0.94) (Fig. S2, Additional file 1). Fitted NOGR values were used for further downstream analysis.

Fig. 1

Auranofin monotherapy. A Image-based quantification of organoid growth rate based on the Viability (V) metric, normalized to timepoint 0 (mean ± SD, n = 2). B Fitted dose–response curves of the mean (n = 2) Normalised Organoid Growth Rate (NOGR) metric for the treatment with Auranofin (500, 909, 1651, 3000 nM) for 120 h. C Normalized Area Over the Curve (AOC) values of the fitted NOGR dose–response curves following Auranofin treatment. D Representative images of organoids treated with vehicle (DMSO) or 909 nM Auranofin for 120 h. Magenta = label-free organoid segmentation by Orbits®; Green = raw cytotox green signal; LU_ = normal pulmonary organoids; NKI_ = lung cancer organoids; PDAC_ = pancreatic ductal adenocarcinoma organoids

AF is selective towards NSCLC and PDAC organoids compared to pulmonary organoids

Organoids were treated with 4 concentrations of AF (500, 909, 1651, 3000 nM) for 120 h and a strong variability in response was observed between the patients (Fig. 1B-C). Notably, the normal pulmonary organoids, specifically LU_46N and LU_51N, demonstrated the highest resistance. They showed no measurable response to AF at concentrations below 1000 nM. Figure 1D presents images of organoids treated with 909 nM AF. These images highlight the precision of the NOGR metric we utilized to categorize the cytostatic effect of the drug at this concentration, where 1 > NOGR > 0. Cytotoxicity could only be achieved in certain tumor organoid lines at higher concentrations of AF, which are more challenging to achieve in the patient as discussed above (> 1500 nM). This highlights the need for combination strategies to enhance the efficacy of AF.

Low CA12 expression is associated with high AF sensitivity

We classified the 9 organoid lines in resistant, intermediate, and sensitive based on the normalized area over the curve (AOC) of the fitted NOGR dose–response curve (Fig. 2A). A higher value indicates a stronger response. Using this classification, we performed a biomarker analysis on the baseline transcriptome of the organoids, which ranked high IGFL1 expression as the strongest biomarker for resistance and low CA12 expression as the strongest biomarker for sensitivity (Fig. 2B). Next, we tested this decision tree on a new PDAC organoid line PDAC_087, with low IGFL2 expression and high CA12 expression (Fig. 2C-D). Based on these markers, we correctly classified PDAC_087 as an intermediate responder, corresponding to its NOGR AOC value below 0.4295 (Fig. 2A). From the partial correlation network for CA12, we identified a positive correlation with NFKB1. Since AF is known to inhibit NF-κB signaling, we performed a correlation analysis of CA12 and NF-κB related genes and found a strong positive correlation with among others the NKFB1, NFKB2 and RELB subunits and a negative correlation with NFKBID and NKRF, which are negative regulators of NF-κB signaling (Fig. 2F). Furthermore, gene sets related to NF-κB signaling were negatively enriched in the sensitive versus resistant organoids, including the NF-κB survival signaling in response to ROS/Hypoxia (Fig. 2G-H). Consequently, cells with impaired NF-κB signaling are more sensitive to AF treatment. Although we are the first to screen AF in NSCLC and PDAC organoids, we and Li et al. have already screened AF in a combined panel of 17 2D cancer cell lines [2, 18]. By combining the AF response classification with publicly available RNA-Seq data from these cell lines [19], we show that all but one resistant cell lines have high expression levels of CA12 (Fig. 3A). ROC analysis further supported the use of CA12 as an accurate biomarker for AF response (sensitive vs. resistant) in these cell lines (AUC: 0.871, p-value: 0.0112) with a sensitivity of 86% and specificity of 100% at a TPM cut-off of 4.645 (Fig. 3B). Based on the comprehensive analysis of organoid and 2D cancer cell line responses to AF, we show that low CA12 expression is a reliable biomarker for predicting sensitivity to AF treatment. To examine the clinical significance of CA12, we analyzed its expression using single-cell RNA sequencing data from publicly accessible databases, focusing on PDAC with 23 patient samples and NSCLC adenocarcinoma with 44 patient samples, specifically within cancer cell populations [11, 12]. In line with our observations in PDOs, we noted a pronounced heterogeneity in CA12 expression among patients, which was also evident at the level of individual patients between cells (Fig. 3C-F). Expression was generally higher in PDAC patients compared to NSCLC patient. Notably, the proportion of CA12-expressing cancer cells within a patient exhibited a strong positive correlation with key components of the NF-κB signaling pathway — NFKB1 (r = 0.73, p < 0.0001), NFKB2 (r = 0.67, p < 0.0005), RELA (r = 0.64, p < 0.011), and RELB (r = 0.51, p < 0.0127) — in PDAC samples (Fig. 3G). This pattern was not replicated in NSCLC samples (Fig. S3, Additional file 1), possibly due to their lower overall expression levels of CA12. Additionally, the quality of the PDAC dataset was markedly superior to that of the NSCLC dataset, which necessitates caution when drawing definitive conclusions from the latter (Fig. S4, Additional file 1).

Fig. 2

Predictive biomarkers for Auranofin response. A Percentile-based classification of Auranofin response into sensitive, intermediate and resistant groups. B Decision tree-based classification of the 3 response groups based on the topmost important features. PDAC_087 was excluded from this analysis to be used as a testing sample. Log2(counts per million reads, CPM) expression values for IGFL1 (C) and CA12 (D). E Partial correlation network for CA12. Grey edges correspond to positive correlation, red edges to negative correlation. The width of the edge is proportional to the absolute partial correlation value of the gene pair. F Positive and negative correlation of CA12 expression with members of the NF-kB signaling pathway. G Gene set enrichment analysis of NF-kB signaling related genesets, following differential expression analysis between the sensitive and resistant organoids. H Up- (red) and down- (blue) regulated genes of the NF-kB survival signaling pathway in the resistant versus sensitive organoids

Fig. 3

CA12 as predictive biomarker in 2D cancer cell lines and expression in patient samples. A CA12 expression values (transcript per million) derived from a publicly available RNA-Seq dataset for 17 2D NSCLC and PDAC cancer cells lines for which the AF treatment response was defined in previous studies. B ROC analysis for the classification of resistant and sensitive (sensitive + intermediate) cell lines to AF treatment, based on CA12 expression values. UMAP overview of cancer cells annotated by patient and CA12 expression for (C) PDAC and D NSCLC adenocarcinoma patient samples. Dotplot representing the fraction of positive cells (%) and mean expression per patient for (E) PDAC and (F) NSCLC patients. G Scatter plot visualising the correlation between the percentage PDAC positive cells for CA12 and NFKB-related genes. The Spearman correlation coefficient and related p-value is plotted. (p < 0.05 indicates significance)

Predictive transcriptome signatures

Besides looking at a single biomarker, we aimed to make predictive transcriptome signatures based on the top correlated genes with the AF NOGR AOC values of the 9 organoid lines. The heatmap in Fig. 4A shows the top 25 negatively (S1) and top 25 positively (S2) correlated genes with the corresponding functional annotation to the Hallmark gene set collection of cluster S1. Interestingly, samples with high p53 pathway activity, IL6_JAK_STAT3 and PI3K_AKT_MTOR signaling appeared to be more resistant to AF, which is in line with previous findings and support the biological relevance of our approach (Fig. 4B). Next, we made two predictive signatures from significantly (p < 0.01) positive correlated genes (n = 86) and negative correlated genes (n = 128) and performed gene set enrichment analysis for PDAC_087 versus the Sensitive, Intermediate and Resistant organoids to determine if we could correctly classify PDAC_087 as an intermediate responder based on our signatures. Compared to the resistant group, PDAC_087 was positively enriched (NES = 2.791, p = 0.003) for genes in the positive correlated gene signature (PCGS) and negatively enriched (NES = -2.324, p = 0.003) for genes in the negative correlated gene set (NCGS) indicating that PDAC_087 will respond to AF (Fig. 4C). In contrast, PDAC_087 was negatively enriched (NES = -1.859, p = 0.003) for genes in the PCGS compared to the sensitive group, and positively enriched (NES = 2.041, p = 0.003) for genes in the NCGS, indicating that it has an intermediate response (Fig. 4D). These results successfully demonstrate the utility of predictive transcriptome signatures in classifying organoids into appropriate AF response categories.

Fig. 4

Predictive signatures for Auranofin response. A Clustered heatmap of the top 25 negative (cluster S1) and top 25 positive (cluster S2) correlated genes with the normalized Area Over the Curve (AOC) values of the fitted dose–response curves for Auranofin in 9 organoid lines, excluding PDAC_087. B Functional annotation of cluster S1 for the Hallmark geneset collection. C-D Volcano plots and normalized enrichment score (NES) of the predictive signatures derived from significantly (p < 0.01) positive (n = 86) and negative correlated genes (n = 128) for PDAC_087 versus the sensitive, intermediate and resistant organoids

Synergy

The objective of our study was to evaluate the effectiveness of AF when used together with eleven different anticancer drugs. For this purpose, we administered each of these drugs in a six-level logarithmic concentration range from 10 to 5000 nM, in combination with AF administered in a four-level concentration range from 500 to 3000 nM. This approach created a synergy matrix of 6 × 4, allowing for an extensive assessment of the combined effects of AF and each anticancer drug. For every combination within this matrix, we quantified the NOGR values. This measurement enabled us to evaluate both the cytostatic (cell growth-inhibiting) and cytotoxic (cell-killing) effects from live-cell imaging as detailed in the materials and methods section and visualized in Figure S1 (Additional file 1). By analyzing these data, we were able to identify the combinations that most effectively induced cell death, pointing to the most promising strategies for combination therapy involving AF and these anticancer drugs. In order to thoroughly determine the degree of combination synergy and select the best model for our study, we compared the commonly-used synergy models HSA, Loewe, Bliss, and ZIP. The results from ZIP and Bliss showed a strong Pearson correlation (r = 0.96), suggesting similar outcomes (Fig. S5A, Additional file 1). However, the conclusions drawn from these models compared to HSA and Loewe varied significantly (Fig. S5B-F, Additional file 1). For instance, the combination of AF with Everolimus was identified as synergistic by the HSA and Loewe models, but only additive by the Bliss and ZIP models in the NKI-120 context. Notably, several combinations did not meet the basic synergy concept (1 + 1 > 2), as illustrated in Fig. S5G-I (Additional file 1). Based on its stringency and accuracy in quantifying true synergism, we selected the ZIP model as the most appropriate for our analysis.

Figure 5A presents the targets for the 11 compounds we selected for our combination screening, and Fig. 5B lists the normalized AOC values of the fitted NOGR curves for each monotherapy across all organoid lines. Notably, AF emerged as the only compound demonstrating selectivity for cancer cells over healthy pulmonary cells. The most pronounced synergistic and cytotoxic effect, as indicated by a high combination sensitivity score (CSS), was observed in combination with the AKT inhibitor MK2206 (Fig. 5C). This effect was particularly significant in the group characterized as intermediate responders to AF and absent in the healthy organoids (Fig. 5D-E). In the highly sensitive organoids, an additive effect was still obtained, resulting in high cell death of the organoids at low concentrations of MK2206 and AF (909 nM) (Fig. 5D), making this a highly potent combination strategy. Interestingly, this was not observed when blocking the PI3K-AKT-MTOR axis upstream (Buparlisib, Fig. S6A, Additional file 1) or downstream (Everolimus, Fig. S6B, Additional file 1) of AKT.

Fig. 5

Auranofin drug combination strategies. A Overview of the 11 therapeutics that were tested in combination with Auranofin. B Heatmap showing the normalized Area Over the Curve (AOC) values of the fitted Normalised Organoid Growth Rate (NOGR) dose–response curves for each therapy and organoid line. C Bubble plot showing the mean ZIP synergy score (bubble size) and combination sensitivity score (CSS, colored heatmap) for each drug and organoid line. D Bubble plot showing the ZIP synergy score (bubble size) and NOGR (colored heatmap) for a concentration range of MK2206 and 909 nM Auranofin for each organoid line. Bubble size: small = ZIP < -10 indicating antagonism; medium = -10 < ZIP < 10 indicating an additive effect; large = ZIP > 10 indicating synergism. NOGR between 1 and 0 indicates a cytostatic response, NOGR < 0 indicates a cytotoxic response. The resistant, intermediate and sensitive groups refer to the Auranofin response classification. E Representative images of organoids (PDAC_044) treated with vehicle (DMSO), 416 nM MK2206, 909 nM Auranofin or the combination for 120 h. Magenta = label-free organoid segmentation by Orbits®; Green = raw cytotox green signal; LU_46N = normal pulmonary organoid; PDAC_044 = pancreatic ductal adenocarcinoma organoid

Investigating the MAPK/ERK pathway, we found that Trametinib (MEK1/2) and ASTX029 (ERK1/2) did not selectively target cancer cells. This indicates that the MAPK/ERK pathway may be essential for the survival of normal epithelial cells, at least in organoids (Fig. 5B). Trametinib showed cytotoxic synergy with 909 nM AF in some tumor organoids at higher concentrations, but this was not selective for cancer cells (Fig. S6C, Additional file 1). Similarly, ASTX-029 demonstrated only a few instances of synergistic interactions (Fig. S6D, Additional file 1), akin to the multi-kinase inhibitor Anlotinib (Fig. S6E, Additional file 1). For IM156 and Palbociclib, a moderate to strong selective synergistic effect was observed when combined with 909nM AF in PDAC_002 and PDAC_060, both intermediate AF responders. However, it's important to note that this synergy was predominantly cytostatic, failing to induce cell death in cancer cells within the nanomolar range (Fig. 6A-B).

Fig. 6

Selected Auranofin drug combination strategies. Bubble plots showing the ZIP synergy score (bubble size) and Normalized Organoid Growth Rate (NOGR, colored heatmap) for a concentration range of A IM156, B Palbociclib, C Cisplatin and D Paclitaxel in combination with 909 nM Auranofin for each organoid line. Bubble size: very small = ZIP < -10 indicating antagonism; small = -10 < ZIP < 10 indicating an additive effect; medium = 10 < ZIP < 20 indicating moderate synergism; large = ZIP > 20 indicating strong synergism. NOGR between 1 and 0 indicates a cytostatic response, NOGR < 0 indicates a cytotoxic response

Finally, a selection of standard of care chemotherapy agents showed moderate to strong cytotoxic synergy when combined with Cisplatin, particularly in the Cisplatin resistant NKI_031 and intermediate responders NKI_120 and NKI_127 (Fig. 6C, Fig. 4B) Notably, this synergistic effect was not observed in the organoid line that exhibited the strongest response to Cisplatin (NKI_142, Fig. 5B). In a similar pattern, synergy with Paclitaxel was observed exclusively in PDAC_060, which was the least responsive to Paclitaxel (Fig. 6D, Fig. 5B). While the concentration range for Gemcitabine was suboptimal due to its strong cytotoxic effect as a monotherapy, a notable observation was made at its lowest concentration. Here, a moderate synergistic effect with 909 nM AF was detected in the most resistant organoids, NKI_031 and PDAC_060, as shown in Figure S6F-G (Additional file 1).