The PDAC-derived cell line BxPC3 (classical subtype) and Panc-1 (basal-like subtype) were cultured in DMEM high glucose with 10% fetal bovine serum and 1% Penicillin-Streptomycin-Glutamine (Sigma-Aldrich, St. Louis, USA) at 37°C, 5% CO₂ in a humidified atmosphere. Mycoplasma contamination was excluded in both cell lines by PCR (MycoScope PCR Detection Kit, Genlantis, San Diego, CA). Cell line authentication was performed by short tandem repeat (STR) profiling using the PowerPlex^® 21 System (Promega, Madison, USA) according to the manufacturer’s instructions.

To generate single cell-derived cell lines (SCDCLs), limited dilution of both parental cell lines was performed in 96-well plates. For Panc-1 two and for BxPC3 four 96-well plates were seeded initially. Each well was examined by phase-contrast microscopy three hours after plating to ensure that only wells harboring a single cell were used for further cultivation. Once the single-cell clones reached approximately 80% confluency in a 96-well plate, they were transferred to a 6-well plate. Upon reaching 80% confulency in the 6-well-plate, the SCDCLs were further transferred to a T-25 flask. After reaching 80% confluency in the T-25 flask, both RNA and protein samples were collected.

To measure the sensitivity to gemcitabine (Sigma-Aldrich, St. Louis, USA), cells of the respective parental cell population growing in log phase (2,000 cells of BxPC3 and 1,500 cells of Panc-1 per well) were seeded in 96-well plates. After a period of 24 hours, the cells were treated with gemcitabine concentrations ranging from 2.5 to 320 nmol/L for 72 hours. After this time, the survival fraction of the cells was determined. For this purpose, cell metabolism was used as a surrogate parameter for viability by using the CellTiter-Blue assay (Promega, Madison, USA). Normalization was based on untreated controls. IC₅₀ and IC_max of parental cell populations were determined using dose-response curves. Subsequently, each SCDCL was treated with gemcitabine in the same experimental set-up at the IC₅₀ and IC_max of the respective parental cell population. Finally, for selected SCDCLs, complete dose-response curves were generated, as performed for the parental cell lines. Each individual measurement was carried out in triplicates and repeated three times.

Cells were seeded at 20,000 cells per well in a 6-well plate(Sarstedt AG & Co. KG, Nümbrecht, Deutschland). Cells were detached from a well by trypsination and counted with a Neubauer chamber (Brand GmbH & Co. KG, Wertheim, Deutschland) every 24 hours. Counting results were normalized to day 1. For each proliferation curve, three independent biological replicates were performed for each time point.

RNA was extracted from each SCDCL using the AllPrep RNeasy Mini Kit (Qiagen N.V., Venlo, Niederlande) as indicated by the manufacturer in the instruction manual. RNA samples were sequenced at Novogene Europe, Cambridge, United Kingdom. A poly-A enrichment and strand-specific library preparation were used. Sequencing was performed on an Illumina Novaseq6000 with S4 flowcell and PE150 length aiming for 30 million reads per sample. The data have been deposited to GEO Accession with the data set identifier GSE232549.

Raw sequencing data (fastq format) was mapped against the human transcriptome (Ensembl GRCh38.103) using kallisto (v0.46.1) and differential expression analysis was performed using sleuth (v0.30.0) (18, 19). Gene set enrichment analysis (GSEA) on b-values (effect sizes estimated by sleuth) was performed using mitch (v1.4.1) against HALLMARK gene sets extracted from the msigdf R package (v7.0) (20).

Protein from each SCDCL was extracted using the EasyPep™ Mini MS Sample Prep Kit (Thermo Fisher Scientific Inc., Waltham, USA). Extracted protein was analyzed using label-free micro-LC tandem mass spectrometry (Ultimate 3000 nHPLC, ThermoFisher & 5600+ Triple TOF, AB Sciex) using data-independent acquisition (DIA). After digestion of non‐labeled protein samples with trypsin, transmitted ions were fragmented and analyzed in the TOF MS Analyzer at high resolution. The raw SWATH data were processed using the software tool DIA-NN v1.7.16 (data-independent acquisition by neural networks) developed by Vadim Demichev et al. (21). The software was used in the high accuracy LC mode with RT-dependent cross-normalization enabled. Mass accuracy, MS1 accuracy, and scan window settings were set to 0, as DIA-NN optimizes these parameters automatically. The ‘match between runs’ function was used to first develop a spectral library using the ‘smart profiling strategy’ from the data-independent acquisition data. The human UniProtKB/swiss-prot database (version 2020/12/6) was used for protein inference from identified peptides (22). Trypsin/P was specified as protease. The precursor ion generation settings were set to peptide length of 7–52 amino acids, the maximum number of missed cleavages to one. The maximum number of variable modifications was set to zero. N-terminal methionine excision and cysteine carbamidomethylation were enabled as fixed modifications. The resulting report file was further processed in the DIA-NN R package for MaxLFQ-based protein quantification (21, 23). A report was generated containing unique proteins (proteins that were not assigned to a group of homologs) that passed the FDR cut-off of 0.01 applied on the precursor level and were identified and quantified using proteotypic peptides only.

The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD042256 (24, 25).

For calculation of the most important proteins that were related to the heterogeneity of SCDCLs we used the R software (v. 4.1.2) along with the packages caret and klaR (26–28). Both packages are available from the CRAN.repository (https://cran.r-project.org/). Using the random forest (RF) algorithm implemented in the caret package the protein data obtained from the MS experiments were used to fit a regression model. Measured response to the IC₅₀ value of gemcitabine of the respective parental tumor population was used as target variable of the model. ‘Backward Feature Elimination’ was utilized to select the most important proteins for the best fit of the regression model. To this end, the ‘Root Mean Square Error’ (RMSE) served as performance measure and the protein set that yielded the best RMSE value was selected from each model run. Using different seed values we took advantage of the random characteristics of the RF algorithm and performed several replicates for both cell lines (10 for Panc-; 30 replicates for BxPC3) of the model runs. In our model, a higher ranking of a protein within a list of the respective run corresponded to a higher relevance for explaining heterogeneous therapy response to IC₅₀ of gemcitabine of the respective parental tumor population. As each run computed an at least slightly different set of ‘optimal’ proteins the results of all model runs were integrated by calculating a total score for each protein as follows: ∑ k = 1 n X k = 1 − a − 1 b (b=“number of proteins in a run”, a=“rank position of the protein in this run”, n=total number of runs). We arbitrarily chose 50 as the upper limit for the number of proteins in the final list because the average number of proteins determined for BxPC3 corresponded to this order of magnitude, whereas the number of proteins determined for Panc-1 was significantly lower.

To measure the sensitivity to JQ1 (APExBIO, Houston, USA), cells of the respective parental cell population or SCDCLs growing in log phase (2,000 cells of BxPC3 and 1,500 cells of Panc-1 per well) were seeded in 96-well plates. After a period of 24 hours, the cells were treated with concentrations ranging from 1.92 to 30,000.00 nmol/L for 72 hours. After this time, the survival fraction of the cells was determined. For this purpose, cell metabolism was used as a surrogate parameter for viability by using the CellTiter-Blue assay (Promega, Madison, USA). Normalization was based on untreated controls. Each individual measurement was carried out in triplicate and repeated three times.

If not stated differently, all analyses were performed using R version 4.1. Responder stratification was performed using stratifyR (v1.0-3) with 2 strata and a fixed total sample size of 0.9 (29). Data handling was performed using the tidyverse package (v2.0.0) including ggplot2 for plotting.

We hypothesized that molecular preconditions of tumor cell subclones within heterogeneous pancreatic cancers correlate with differential response to therapy and could be targeted to modify treatment response ( Figure 1A ). Hence, we established single cell-derived cell lines (SCDCLs) of the parental cell populations from the classical differentiated cell line BxPC3 and basal-like cell line Panc-1 by limiting dilution. Twelve SCDCLs of BxPC3 and 14 SCDCLs of Panc-1 were generated as schematically shown in Figure 1B . The time period after single cell sorting until 80% confluency in a 6-well culture plate was considerably different between individual SCDCLs of both parental cell lines. Among the SCDCLs of BxPC3, the first one reached confluency after 32 days and the last one after 62 days ( Figure 1C ). The time to 80% confluency of the SCDCLs of Panc-1 ranged from 30 to 48 days ( Figure 1D ).

Cell morphology and growth patterns of the growing colonies differed between distinct SCDCLs of both cell lines as exemplified by the phase contrast images shown in Figures 1E, F . Although the BxPC3 SCDCLs all tended to grow in rather dense formations, we also observed a more elongated shape of the cells in some of the SCDCLs, while others grew much more cuboidal ( Figure 1E ). Notably, a higher number of spindle-shaped cells were observed in the SCDCLs of Panc-1 in some clonal cultures, indicating a more mesenchymal phenotype. Other SCDCLs of Panc-1, however, grew in more cobblestone-like formations and resemble more of a cuboid shape, which indicates a more epithelial phenotype, confirming earlier observations ( Figure 1F ) (17).

To test our central hypothesis, that different SCDCLs of the same parental tumor cell population of PDACs bear distinct intrinsic molecular profiles that determine the individual response to chemotherapy, each SCDCL was tested for its individual response to gemcitabine in vitro. First, we determined the half-maximal inhibitory concentration (IC₅₀) and a concentration close to the maximal inhibitory concentration (IC_max) of the BxPC3 and Panc-1 parental cell populations. Dose-response curves to gemcitabine were generated for both parental cell lines in the concentration range from 2.5 to 320 nmol/L. The parental Panc-1 population (IC₅₀ = 43 nmol/l) proved to be much more resistant to gemcitabine than parental BxPC3 population (IC₅₀ = 9.6 nmol/L) as measured by their respective IC₅₀ ( Supplementary Information S1 ).

Subsequently, each SCDCL of both parental cell populations was treated with the half-maximal inhibitory concentration (IC₅₀) and a concentration close to the maximal inhibitory concentration (IC_max) of the respective parental population. We observed a highly variable response of distinct SCDCLs of both cell lines BxPC3 and Panc-1 ( Figures 2A, B ; Supplementary Information S2 ). Heterogeneity in response was substantially higher among the BxPC3 SCDCLs compared to the Panc-1 SCDCLs.

In detail, the SCDCLs of BxPC3 showed a highly variable response. The most resistant SCDCL B2D9 had a survival rate of 0.67, while this was only 0.29 for the most sensitive SCDCL B3C10 when treated with IC₅₀ of the parental cell population (p<0.009). There was a continuum of response rates to gemcitabine treatment between these extremes, although one might suspect a greater increase from the fourth most resistant SCDCL onwards ( Figure 2A ).

Among the SCDCLs of Panc-1, there was also a heterogeneous response, although the differences between the most resistant SCDCL P4E2 with a survival rate of 0.8 and the most sensitive SCDCL P4B9 with 0.66 (p<0.006), when treated with IC₅₀ of the parental cell population, were considerably smaller. The SCDCLs with survival rates between these extremes formed a much denser continuum compared to BxPC3.

From the most sensitive and the most resistant SCDCL of both parental cell lines, complete dose-response curves were subsequently established, corresponding to the experimental design described above. These curves showed significant differences in the resistance profile between BxPC3 SCDCls B2D9 and B3C10 in concentration ranges from 5 nM to 80 nM (p<0.05) ( Figure 2C ). Among the Panc-1 SCDCLs P4E2 and P4B9, a significant difference was only observed at concentration ranges from 2.5 nM to 5nM (p<0.05) ( Figure 2D ).

Next, we also tested for the other SCDCLs whether the resistance to gemcitabine occurred randomly at the IC₅₀ of the respective parental population, or rather indicated a more resistant behavior of the SCDCL in general. The response rates of the SCDCLs at the IC₅₀ and IC_max tended to correlate for both cell lines (SCDCLs of BxPC3: r=0.4828; p<0.001; SCDCLs of Panc-1 r=0.3089; p<0.001) ( Supplementary Information S3 ). Thus, we concluded that a higher resistance to gemcitabine at the IC₅₀ of the respective parental population tends to reflect a higher resistance to gemcitabine of the respective SCDCL in general.

Next, we tested whether there is a correlation between the time to confluency after single-cell sorting, which could be an indicator of better adaptive behavior and treatment response. For the SCDCLs of both cell lines, no clear association between the time to confluency and response to treatment was observed ( Supplementary Information S4 ).

Gemcitabine, a deoxycytidine analog, causes inhibition of DNA chain elongation in addition to several other processes (30). To explore a potential correlation between proliferation rate and gemcitabine sensitivity, we determined population doubling times (PDT) of the most resistant SCDCLs, most sensitive SCDCLs and the parental tumor population of both cell lines. The most resistant SCDCL of BxPC3 (B2D9: PDT= 44.27h; CI95%: 41.29h to 47.72h) had a higher PDT in comparison to the most sensitive SCDCL (B3C10: PDT= 29.28h; CI95%: 26.99h to 32h) and the parental BxPC3 population (PDT: 31.47h; CI95%: 29.05h to 34.31h). For further validation, we additionally determined the PDT of the second and third most gemcitabine resistant SCDCL. However, PDTs of these SCDCLs were similar compared to the most sensitive SCDCL B3C10 and the parental cell line BxPC3 (B2F8; PDT= 34.13h; CI95%: 30.72h to 38.38h; B1G3: PDT= 32.82h; CI95%: 29.13h to 37.59h) ( Figure 2E ). There were no significant differences between the parental Panc-1 population (PDT= 26.86h CI95%: 25.18h to 28.79h), the most sensitive (P4B9: PDT= 29.91h; CI95%: 27.31h to 33.06h) and most resistant SCDCL (P4E2: PDT= 29.42h; CI95%: 26.33h to 33.34h) ( Figure 2F ).

We concluded that the observed differences in response to gemcitabine are a reflection of molecular idiosyncrasies of the individual SCDCLs that are independent of their intrinsic proliferation rates. Thus, we comprehensively profiled the transcriptome and proteome of individual SCDCLs of both parental cell populations to elucidate the molecular basis of the response of SCDCLs to gemcitabine.

To determine whether the distinct treatment phenotypes of the SCDCLs are related to transcriptomic heterogeneity, we performed mRNA sequencing of all 12 SCDCLs of BxPC3 and all 14 SCDCLs of Panc-1.

Unsupervised Principle Component Analysis (PCA) on the 1,000 most variably expressed genes showed a clear clustering of the SCDCLs according to their respective sensitivity to gemcitabine treatment for BxPC3, while this clustering was less clear for Panc-1 ( Figures 3A, C ; Supplementary Information 5A, C ).

For data evaluation, SCDCLs of each parental cell population were divided into two groups according to their respective treatment response, i.e., a more resistant group (resG) and a more sensitive group (sensG). The cut-off value of the survival rate for the assignment of each SCDCL to the respective group was calculated for treatment at the respective IC₅₀ as well as IC_max according to Reddy et al. (29). For the SCDCLs of BxPC3, the cut-off values based on the IC₅₀ and IC_max were 0.44 and 0.18, respectively. The cut-off values for Panc-1 were 0.73 and 0.49 based on the IC₅₀ and IC_max, respectively.

Next, we aimed to identify differentially expressed genes between SCDCLs of the resG compared to the sensG. When stratified according to the IC₅₀ value, the SCDCLs of BxPC3 showed 159 differentially expressed genes (q < 0.1; 106 up-regulated in resG, 53 up-regulated in sensG) ( Supplementary Information 5B ). When stratified by IC_max, 753 genes were differentially expressed (q < 0.1; 327 up-regulated in resG, 426 up-regulated in sensG) ( Figure 3B ). Consistent with less heterogeneity in response to gemcitabine, Panc-1 SCDCLs showed fewer differentially expressed genes. When stratified by IC₅₀ value 98 genes were differentially expressed (q < 0.1; 37 up-regulated in resG, 61 up-regulated in sensG) ( Supplementary Information 5D ), and when stratified by IC_max 149 genes were differentially expressed (q < 0.1; 74 up-regulated in resG, 75 up-regulated in sensG) ( Figure 3D ). Measured by the number of differentially expressed genes in the SCDCLs of the respective parental cell populations, we conclude that there is less transcriptional heterogeneity between SCDCLs of the basal-like cell line Panc-1 compared to the SCDCLs of the classical cell line BxPC3. Strikingly, the lower transcriptional heterogeneity is reflected in the lower heterogeneity of response to gemcitabine treatment.

Next, gene set enrichment analysis was performed to identify the differentially regulated signaling pathways between the resG and the sensG ( Figures 4A, C ).

The SCDCLs of both parental cell populations showed a continuum between the most resistant and most sensitive SCDCL in their response to gemcitabine ( Figures 2A, B ). Hence, we hypothesized a rather continuous differential regulation of relevant signaling pathways rather than on-off effects. We performed additional GSEA between the three most resistant SCDCLs (BxPC3: B2D9, B2F8, B1G3; Panc-1: P4E2, P3D10, P1C3) and the three most sensitive SCDCLs (BxPC3: B3C10, B2G7, B3F8; Panc-1: P4B9, P4B5, P3D2) of both, BxPC3 and Panc-1, to be able to generate a better discriminatory power of differentially regulated pathways ( Figures 4B, D ).

The pathways MYC targets V1, MYC targets V2 as well as E2F targets from the Molecular Signatures Database (MSigDB) were enriched in the sensitive SCDCL of BxPC3 across all comparisons. Interestingly, exactly these pathways were enriched in Panc-1 in the group of the resistant SCDCLs across all comparisons. In contrast, enrichment of EMT correlated with resistance to gemcitabine in BxPC3 SCDCLs across all comparisons, while Panc-1 SCDCLs showed no resistance-associated enrichment of this pathway.

As described above, we hypothesize that relevant signaling pathways are gradually differentially regulated in the SCDCLs. Hence, particular attention was paid to pathways that did not show significant differential expression in the resG and sensG comparisons, but in contrast, were differentially regulated in the comparison of the three most resistant versus the three most sensitive SCDCLs. Since a positive correlation was observed between the response rates of the SCDCLs to therapy at the IC₅₀ ( Figures 4A–D ) as well as the response rate at the IC_max ( Supplementary Information S6A–D ), the differentially enriched pathways overlapping in the different comparisons were further considered. Overlaps of enriched gene sets between all comparisons are shown in a tabular overview in Supplementary Information S6E, F . Consequently, in BxPC3 SCDCLs, enriched signaling pathways EMT, TNF signaling via NfKB, and IL2STAT5 signaling correlated with more resistant behavior to gemcitabine. In addition, several other inflammatory signaling pathways appear to be associated with resistance. An overview of the differentially regulated pathways according to IC_max classification is shown in Supplementary Information 6 . In contrast to BxPC3, no additional enriched pathways were found in the SCDCLs of Panc-1, when comparing the three most resistant and the three most sensitive SCDCLs. This result is in line with the lower transcriptional heterogeneity among SCDCLs of Panc-1.

For a comprehensive understanding of the biological processes associated with the heterogeneity of SCDCLs and their distinct intrinsic resistance profiles to gemcitabine, we next sought to profile their individual proteomes. As commonly known, mRNA expression levels do not necessarily reflect the respective protein expression levels (31). To obtain a more comprehensive picture, we therefore aimed to identify proteins that were associated with the heterogeneous response to gemcitabine of individual SCDCLs by mass spectrometry analyses. Using the feature selection with the random forest approach, we extracted protein signatures that are associated with the response to gemcitabine of each individual SCDCL. Extracted proteins were ranked according to their predictive value for treatment response. This approach reflects the fact that different protein compositions may be similarly important for adaptation to the IC₅₀ target variable

For BxPC3, we extracted 50 proteins that were associated with the functional heterogeneity of response to gemcitabine of individual SCDCLs ( Supplementary Information S7 ). Of these, overexpression of 21 proteins was associated with a poorer response to gemcitabine, whereas overexpression of 29 proteins was associated with a better response to gemcitabine ( Figure 5A ). TNF receptor superfamily member 6b (TNFRSF6B) was ranked highest among the proteins extracted for BxPC3 and associated with a poorer response. In addition, DNA activity-influencing proteins such as bromodomain 3 (BRD3) and high mobility group nucleosome binding domain 5 (HMGN5) were highly ranked in the identified set of proteins associated with a poorer response. Among others, HMGN5 was also expressed significantly stronger at the mRNA level in the resG compared to the sensG (p=0.046).

For Panc-1, we again identified substantially fewer proteins whose expression level was associated with response to treatment of individual SCDCLs which is in line with the lower transcriptional heterogeneity described above. A set of 18 proteins was identified that was associated with the heterogeneity of response to gemcitabine of individual SCDCLs ( Supplementary Information S7 ). Of these, overexpression of six proteins was associated with a poorer response to gemcitabine, whereas overexpression of 12 proteins was associated with a better response to gemcitabine ( Figure 5B ). Autocrine Motility Factor Receptor (AMFR) was ranked highest among the proteins extracted for Panc-1 ( Supplementary Information S7 ).

JQ1 is an inhibitor of the bromodomain and extraterminal family of proteins (BET) with the highest selectivity for BRD4 (32).

In our pathway enrichment analyses based on transcriptomics described above, the more gemcitabine-resistant SCDCLs of the classical cell line BxPC3 showed enrichment of EMT, TNF signaling via NfKB, and IL2STAT5 signaling. The proteome analyses of the same SCDCLs identified two pivotal proteins in the extracted protein signature, i.e., (i) BRD3 which is a family member of the BET proteins and (ii) the TNF receptor TNFRSF6B whose gene possess a super-enhancer in multiple myeloma cells (33).

Pathway enrichment analyses based on transcriptomics of the gemcitabine-resistant SCDCLs of the basal-like cell line Panc-1 showed enrichment of MYC signaling, i.e., gene sets MYC targets V1 and MYC targets V2.

JQ1 suppresses cell proliferation through several signaling pathways such as TNFA_Signaling_via_nfkb, L2_STAT5_SIGNALING, MYC signaling as well as multiple inflammatory transcriptional programs in pancreatic cancer (33–37).

As described above, we intended to identify molecular preconditions of tumor cell subclones that could be targeted to overcome treatment resistance of current clinical standard therapy ( Figure 1A ). In our model system, these subclones are reflected by SCDCLs of heterogeneous parental cell populations which show differential response to treatment with gemcitabine. Based on our molecular findings, we hypothesize that gemcitabine-resistant SCDCLs could be specifically targeted by inhibition of proteins of the BET family in both, the classical cell line BxPC3 and the basal-like cell line Panc-1. Therefore, we tested the specific anti-proliferative effect of JQ1, which is an effective inhibitor of the BET proteins, in our SCDCL model system (32).

In general, the parental population of Panc-1 (IC₅₀: 679 nM CI95%: 496.5 – 929.1) was more resistant to JQ1 treatment compared to the parental BxPC3 population (IC₅₀: 184.3 nM; CI95%: 146.3 - 232.2). Next, we generated dose-response curves for JQ1 of the parental cell population, the most gemctabine-resistant SCDCL and the most gemcitabine-sensitive SCDCL for both cell lines, BxPC3 and Panc-1.

In BxPC3 the most gemcitabine-resistant SCDCL B2D9 was the substantially more sensitive to treatment with JQ1 (IC₅₀: 48.36 nM; CI95%: 27.68 - 84.48). compared to the most gemcitabine-sensitive SCDCL B3C10 (IC₅₀: 95.63 nM CI95%: 73.16 - 125) ( Figure 5C ). Differential response was significant in concentration ranges from 9.6 nM to 48 nM (p<0.02). The parental BxPC3 population was significantly more resistant to JQ1 treatment (IC₅₀: 184.3 nM; CI95%: 146.3 - 232.2) compared to both derived SCDCLs.

The most gemcitabine-resistant SCDCL of Panc-1 again was most sensitive when treated with JQ1 (P4E2: IC₅₀: 471.8 nM CI95%: 338.4 – 658) in comparison to the most gemcitabine-sensitive SCDCL P4B9 (IC₅₀: 1590 nM CI95%: 667 – 3791). The differential response was significant in concentration ranges from 240 nM to 6000 nM (p<0.02). The parental population Panc-1 was in-between these two SCDCLs (IC₅₀: 679 nM CI95%: 496.5 – 929.1) ( Figure 5D ).

In conclusion, we showed that gemcitabine-resistant subclones, i.e., SCDCLs, of the heterogeneous parental cell populations of both the classical cell line BxPC3 and the basal-like cell line Panc-1 can be specifically targeted using the BET inhibitor JQ1.