Single-Cell–Derived Primary Rectal Carcinoma Cell Lines Reflect Intratumor Heterogeneity Associated with Treatment Response

Colorectal cancer is the third most commonly diagnosed cancer in males and the second in females worldwide (1). About 30% of these carcinomas are located in the rectum (2). Preoperative chemoradiotherapy (CRT) followed by radical surgical resection represents the current standard of care for locally advanced rectal cancer (UICC II and III; ref. 3). The response of individual tumors to preoperative CRT, however, is extremely diverse. It ranges from complete response without any residual viable tumor cells detectable on histopathologic examination of the surgical specimen to the complete absence of local tumor regression. Approximately 25% of patients show primary resistance, that is, no clinical response (4). This poses a clinical dilemma, because patients with a priori–resistant tumors could either be spared exposure to CRT and its severe side effects, treated by alternative, intensified protocols, or selected to undergo immediate surgery. On the other hand, patients for which a complete response could be definitively predicted might not need extensive surgery (5) and could be spared from surgery-associated morbidity and mortality (6). However, identifying reliable response predictors from bulk pretreatment tumor biopsies has proven to be very challenging (7).

There has been emerging evidence over the past few years that intratumor heterogeneity (ITH) is a key determinant of tumor biology, treatment response, and ultimately patient survival (8). Genomic instability is one of the hallmarks of tumorigenesis and the basis of ITH (9). It results in the continuous acquisition of genetic aberrations in tumor cell populations, the generation of subclones with distinct aberration profiles, and the activation of specific signaling pathways; this results in a disparate susceptibility to therapeutic interventions (8).

Comprehensive deconvolution of the clonal composition of primary tumors can be achieved by single-cell sequencing (10). In fact, single-cell RNA sequencing (RNA-seq) of colorectal cancers revealed remarkable transcriptional heterogeneity previously unappreciated by bulk tumor sequencing (11). Recent studies using multiregion DNA sequencing revealed that ITH in colorectal cancer is likely to be generated by branching evolutionary processes of tumor clones (12).

We hypothesize that the genetic diversity of rectal carcinomas, that is, ITH, is a major determinant for treatment resistance and resistant subclones. Subclones that acquire resistance under treatment and clonally expand may already be present in treatment-naïve tumors. To explore the impact of ITH and the presence of distinct subclones on treatment response, we generated an in vitro model for clonal analysis of both genetic composition and response. We established single-cell–derived cell lines (SCDCL) from a treatment-naïve biopsy of a rectal carcinoma from a patient with poor clinical response to subsequent CRT.

Tissue was collected after approval by the ethics committee of the University of Connecticut (Waterbury, CT; IRB Number: IE-10-240-3) and xenografted in athymic NCR nu/nu mice (Charles River Laboratories). For details see Supplementary Data S1.

After explantation, the xenograft was dissociated and cell lines were established as described in detail in Supplementary Data S2. Cells were checked for Mycoplasma contamination by PCR (MycoScope PCR Detection Kit, Genlantis). Cell line authentication was verified by short-tandem repeat (STR) profiling using the AmpFLSTR Identifiler PCR Amplification Kit (Thermo Fisher Scientific) according to the manufacturer's instructions.

To measure the sensitivity to γ radiation, cells growing in log phase were seeded as single cells in 6-well plates and irradiated 24 hours after seeding with a single dose of 4 Gy of γ radiation using a ¹³⁷Cs irradiator (Mark 1 irradiator, JL Shepherd & Associates). For CRT, cells growing in log phase were seeded as single-cell suspensions in 6-well plates and 5-FU (Sigma-Aldrich) was added to a final concentration of 3 μmol/L 8 hours after cell seeding. Sixteen hours later, cells were irradiated with a single dose of 2 to 8 Gy of γ radiation. Treatment with 5-FU was stopped 8 hours after irradiation by a medium exchange and cell survival was evaluated by a standard colony formation assay and subsequent calculation of the respective survival fraction (SF). Nonirradiated cultures with or without treatment with 5-FU were used for normalization. Cells were grown for 8 days after treatment, fixed with 70% ethanol, dried, and stained with Hemalaun (Merck). Colonies containing >50 cells were scored as survivors. Each experiment was performed in triplicate and repeated three times. Statistical analyses were performed with an unpaired Student t test for irradiation experiments with 4 Gy and two-way ANOVA for dose–response experiments using GraphPad Prism.

To measure the sensitivity to 5-FU alone, 1,500 cells growing in log phase were seeded per well cells in 96-well plates and 5-FU (Sigma-Aldrich) was added 24 hours after seeding to a final concentration in a range of 0–150 μmol/L. After 72 hours, cell metabolism was measured by CellTiter-Blue (Promega). Each experiment was performed in triplicates and repeated three times.

Cells were transiently transfected with 100 ng of the appropriate dual-luciferase assay plasmid from the TCF/LEF Reporter Assay Kit (Qiagen) using Lipofectamine 3000 (Invitrogen) delivered in Opti-MEM Reduced Serum Medium (Gibco) according to the provided TCF/LEF Reporter Assay Protocol (Qiagen). Cells were transfected for 72 hours. Cell lysates were prepared for measurement using the Dual-Luciferase Reporter Assay System (Promega), following the provided instructions. Measurements were taken with a Centro XS³ LB 960 Microplate Luminometer (Berthold Technologies). Fold change was calculated according to a technical note from Promega.

IHC was performed using the following primary antibodies: monoclonal anti-cytokeratin 20 antibody, clone Ks20.8 (1:200, #M7019, Dako) and monoclonal anti-vimentin antibody, clone SP20 (1:400, #RM-9120-s, Thermo Fisher Scientific). For antigen retrieval, sections were incubated in a water bath at 95°C in either Target Retrieval Solution, pH 9.0 (Tris/EDTA, 1:10; #S2367, Dako) for 20 minutes (cytokeratin 20) or Target Retrieval Solution, pH 6.1 (citrate, 1:10; S1699, Dako) for 40 minutes (vimentin). Detection was done using the EnVision Detection System, Peroxidase/DAB, Rabbit/Mouse (# K5007, Dako) according to the manufacturer's instructions. For Periodic Acid Schiff (PAS) reaction tissue sections were deparaffinized, rehydrated, and incubated in periodic acid (1% w/v, freshly prepared; # 10450-60-9, Merck) for 5 minutes followed by incubation in Schiff's reagent (catalog no. 3952016, Sigma) for 5 minutes. Counterstaining was done with hematoxylin.

Cells were grown on coverslips, washed twice with PBS, and fixed with ice cold methanol for 15 minutes at −20°C. Cells were washed with PBS, blocked with 5% goat serum and 0.3% Triton X-100 in PBS for 60 minutes at room temperature, and incubated with the primary antibody (1:400, cytokeratin 20, #13063, Cell Signaling Technology) in PBS supplemented with 1% BSA and 0.3% Triton X-100 at 4°C overnight. Cells were washed with PBS and incubated for 2 hours with the fluorochrome-conjugated secondary antibody (1:500, Alexa Flour 568 goat anti-rabbit, #A1011, Invitrogen). After washing with PBS, coverslips were mounted with ProLong Gold antifade reagent with DAPI (Life Technologies).

Preparation of metaphase chromosome suspension, slide pretreatment, slide denaturation, hybridization, detection, and imaging were conducted as described previously (13).

Bacterial artificial chromosome (BAC) contig assembly, DNA extraction, nick translation, and precipitation as well as sequential hybridization, detection, imaging, and signal counting were performed as reported previously (14). A total of 500 nuclei were analyzed for the parental cell population and 651 for the tumor biopsy. A total of 300 nuclei were analyzed for each SCDCL.

DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen) as instructed by the manufacturer. One μg sample DNA and sex-matched reference DNA was labeled using the SureTag DNA Labeling Kit (Agilent) according to the manufacturer's instructions (protocol version 7.3). Labeled DNA was subsequently hybridized on SurePrint G3 human 4 × 180K CGH Arrays (Agilent) for 24 hours. Slides were processed on an Agilent G2565BA scanner. Agilent Technologies' Feature Extraction (version 10.7.3.1) was used for data quality control and extraction. Data were visualized and analyzed using Nexus 7.5.

RNA was extracted from each cell line using the RNeasy Mini Kit and on-column DNAse Treatment (Qiagen) according to the manufacturer's instructions. RNA integrity was checked using the 2100 Bioanalyzer (Agilent). Sequencing libraries targeting mature mRNA and lncRNA were generated using the TruSeq RNA Kit (Illumina). Samples were sequenced on a single lane of the Illumina HiSeq 2500, with a target depth of ≥25 million reads per sample.

RNA-seq analysis was conducted using the CCBR RNASeq pipeline (https://github.com/CCBR/Pipeliner). For RNA-seq gene expression analysis, contaminating adapter sequences were removed using Trimmomatic v0.36 (15). Trimmed reads were mapped to the GRCh37 human reference genome using STAR v2.5.3 (16) run in 2-pass mode. RSEM v1.3.0 (17) was then used for gene-level expression quantification, and DESeq2 v1.18.1 (18) was used to assess differential expression with FDR correction using the Benjamin–Hochberg method for multiple tests. Only genes with a counts per million (CPM) of at least 0.5 in a minimum of 3 samples were carried forward for differential expression analysis.

Read- and alignment-level quality was assessed using MultiQC v0.9 (http://multiqc.info/) to aggregate QC metrics from FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), FastQ Screen (https://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/), Picard, RSeQC (http://rseqc.sourceforge.net/), and Trimmomatic.

Pathway enrichment analysis of the top 1,000 genes based on P value for each pairwise differential comparison was performed using ToppGene Suite (http://toppgene.cchmc.org; ref. 19) and the Ingenuity Pathway Analysis Tool (IPA; Qiagen). Gene set enrichment analysis (GSEA) was performed using preranked GSEA (20) with default settings and the MsigDB (20) canonical pathway (CP) gene sets. For this analysis, genes were ranked using a t-statistic to produce a list sorted from the most significant upregulated gene down to the most significant downregulated gene. For each preranked GSEA analysis, we used the Enrichment Map (21) plugin in Cytoscape v3.6.1 (22) to collapse pathways based on shared gene membership and functional overlap. All pathway and gene set enrichment P values were FDR-corrected.

DNA was extracted (i) from the tumor biopsy using the DNeasy Blood and Tissue Kit (Qiagen), (ii) from the FFPE-fixed adjacent tissue as reported elsewhere (23), and (iii) from the cell lines using the Gentra Puregene Kit (Qiagen) according to the manufacturer's instructions. cDNA libraries were prepared using the Agilent SureSelectXT Human All Exon V5 plus UTR target enrichment kit. These samples were pooled and sequenced on an Illumina HiSeq2500 with TruSeq V4 chemistry.

Alignment and Variant Calling was performed following the CCBR pipeliner workflow (https://github.com/CCBR/Pipeliner). Reads were trimmed using Trimmomatic v0.36 (15) and mapped to the hs37d5 version of the human reference genome (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz) using BWA-mem v07.15 (24). BAM files were processed using Samtools v1.3 (http://www.htslib.org/), and duplicates were marked using Picard v2.1.1 (http://broadinstitute.github.io/picard/). GATK v3.5.0 (25) was used to perform indel realignment and base recalibration. Read- and alignment-level quality control visualization was performed using MultiQC v1.1 (http://multiqc.info/) to aggregate QC metrics from FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), FastQ Screen (https://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/), Picard, bamtools stats (http://github.org/pezmaster31/bamtools), and Trimmomatic (15). Variant calling was performed in paired mode with MuTect2 (26), and variants were annotated using Oncotator v1.9.1.0 (http://portals.broadinstitute.org/oncotator/).

Copy number was estimated from the exome sequencing data using Sequenza v2.1.2 (27) and FREEC v11.5 (28) jointly. First, Sequenza was run to estimate ploidy and purity, and then these ploidy and purity estimates were used as fixed values in the FREEC analysis to call absolute copy number of all genomic segments. Ploidy and purity estimates from Sequenza were then used as priors to estimate cellular fraction and subclone clustering for each SNP and INDEL using Pyclone v0.13.1 (29). Clusters containing <3 mutations were excluded from further analysis and discussion because they had wide 95% confidence intervals (CI) and/or contained no nonsilent mutations.

We collected treatment-naïve tumor biopsies and adjacent normal tissue prior to CRT from a 59-year-old patient with an adenocarcinoma of the lower rectum (uT3, uN1, M0); the patient responded partially to preoperative CRT (regression grade 2 according to Becker; ref. 30) but rapidly developed liver metastases and local recurrence. The clinical history is summarized in Fig. 1A and Supplementary Data S3. To establish an in vitro model that recapitulates ITH in the primary tumor, we aimed to propagate SCDCLs. The experimental workflow is summarized in Fig. 1B. First, we generated a patient-derived xenograft (PDX). After subcutaneous implantation of a tumor biopsy into an athymic nude mouse, substantial tumor growth was observed after 25 days. The tumor was explanted 16 weeks after initial engraftment at a size of 544 mm³ and examined by hematoxylin and eosin (H&E) staining, IHC, and Periodic Acid Schiff (PAS) reaction (Fig. 1C–E). The primary tumor presented as a moderately differentiated colorectal adenocarcinoma with some gland formation and focal PAS-positive mucin production, while the PDX showed a solid growth pattern without gland formation. By IHC, both the primary tumor and PDX were positive for CK20 as a characteristic marker for colorectal epithelial cells and negative for the stromal cell marker vimentin. The PDX was subsequently reimplanted into a second athymic nude mouse to enable further growth.

The secondary PDX was explanted after having reached 608 mm³ at 7 weeks and converted to a cell line employing the Rho-associated kinase (ROCK) inhibitor Y-27632, as described for the propagation of epithelial cells by Schlegel and colleagues (31). Contaminating murine cells were depleted. The patient-derived rectal cancer cell line, termed WHrec1, showed a typical epithelial morphology (Fig. 1F) and, in line with the PDX, expressed moderate to high levels of CK20 (Fig. 1G). The STR profile is provided in Supplementary Data S4.

We analyzed the karyotype of the parental cell population by spectral karyotyping (SKY) and array-based comparative genomic hybridization (aCGH), which revealed a near-triploid karyotype with numerous chromosomal imbalances typical for colorectal cancer (Fig. 1H). aCGH confirmed the loss of chromosomes 1p, 3p, 6, 8, 10, 15, 17p, 18, 21, and Y, and gains of 1q, 19p, 20, and X (Supplementary Fig. S1A).

We subsequently performed single-cell cloning by limiting dilution and established 23 SCDCLs from the parental cell population. The flow chart in Supplementary Fig. S1B summarizes the analytical procedures. The SCDCLs revealed distinct growth patterns and kinetics (Supplementary Fig. S2A and S2B). After single-cell cloning, the time until each SCDCL reached confluency in a T25 tissue culture flask ranged from 33 to 64 days.

We tested the central hypothesis that individual clones of the primary tumor respond differently to CRT. To accomplish this, we selected 16 SCDCLs that represent the spectrum of disparate growth kinetics after single-cell cloning (Supplementary Fig. S2B) and the parental cell population and exposed them each to 4 Gy of γ radiation. The survival fraction (SF) was measured by colony formation assays. We observed a profoundly different sensitivity to radiation across the various SCDCLs (Fig. 2A). The most sensitive SCDCL (2F4) had a SF4 of 0.07, whereas the most resistant SCDCL (1G10) had a SF4 of 0.18 (P = 0.0012). The parental cell population was between these two outermost SF4 values (0.12). Next, we established dose–response curves of the parental cell population as well as the most resistant and most sensitive SCDCL, which revealed significant differences of the SF between SCDCLs 1G10 and 2F4 over a dose range of 2 to 8 Gy (P < 0.0001; Fig. 2B). The addition of 3 μmol/L 5-FU revealed comparable results (P = 0.0003; Fig. 2C). Treatment with 5-FU alone did not reveal significant differences between the SCDCLs 1G10 (radioresistant; IC₅₀: 7.4 μmol/L; 95% CI: 6.9–8.1) and 2F4 (radiosensitive; IC₅₀: 7.9 μmol/L; 95% CI: 7.2–8.3) in the IC₅₀ value (Fig. 2D).

To determine whether the proliferation rates of the SCDCLs influence response to radiation, we measured their respective growth rates. No differences in growth rates were found between SCDCLs 1G10 and 2F4 and the parental cell population (Fig. 2E). There was also no statistically significant correlation between growth kinetics after single-cell cloning and the SF4 (Supplementary Fig. S2C). We conclude that the observed differences in response are therefore a reflection of the idiosyncrasies of the individual SCDCLs that are independent of their intrinsic proliferation rates. Assuming that DNA repair capacity might be involved in the profoundly different sensitivity of the SCDCLs to irradiation, we used the alkaline Comet Assay to quantify the level of DNA damage repair. Even though we observed a trend toward slower DNA damage repair in the sensitive SCDCL (2F4) compared with the resistant SCDCL (1G10), the differences did not reach statistical significance (Supplementary Fig. S2D). Hence, we used comprehensive genome and transcriptome profiling to elucidate the molecular basis of response of SCDCLs to CRT.

To determine whether the distinct treatment phenotypes of the SCDCLs are related to genetic heterogeneity, we performed whole-exome sequencing (WES; Supplementary Data S5). WES showed TP53 (p.G245S) and KRAS (p.G12S) missense mutations and a 32 bp deletion in CTNNB1 (p.VSHWQQQSYLDS22del) in the tumor biopsy, the parental cell population, and both SCDCLs at a cellular fraction of 1.0. No mutations were detected in BRAF and APC. The large in-frame CTNNB1 deletion occurred in a hotspot region (Supplementary Fig. S2E), that has been described in various tumor entities (32).

To infer clonal population structure from WES data, we used PyClone (29). PyClone analysis identified a total of 12 clusters characterized by distinct sets of mutations; three of the clusters had fewer than three mutations and were excluded from further analysis (Fig. 3A; Supplementary Data S6).The remaining 9 clusters were present at different frequencies across the samples (Fig. 3A). Cluster 0 was detected at high cellular fractions across all samples, and contained the clonal TP53, KRAS, and CTNNB1 mutations. Of particular interest is cluster 7, which is characterized by mutations in ZWINT and SV2A. This cluster was detected at very low cellular fractions in the tumor biopsy and the sensitive SCDCL (2% and 0.5%, respectively), but present at intermediate cellular fraction (46%) in the parental cell population and an elevated cellular fraction (63%) in the resistant SCDCL. This suggests a role of this cluster and its associated mutations in treatment resistance. The complete results for the PyClone analysis, including all mutations and inferred clusters, are summarized in Supplementary Data S6.

To identify mutations associated with treatment response, we next looked for differences in mutation prevalence in the resistant versus sensitive SCDCLs by comparing cellular fractions of mutations in the resistant versus sensitive SCDCLs along with the tumor biopsy and parental cell population. Mutations in 16 genes were found at a higher cellular fraction in the resistant SCDCL (1G10), but at reduced cellular fraction in the sensitive SCDCL (2F4), relative to the parental cell population, respectively. Again, this analysis revealed that mutations in SV2A (p.E675A) and ZWINT (p.K74E), which define the resistant enriched PyClone cluster 7, were detected only in the parental cell population and the resistant SCDCL. In fact, the cellular fractions of both mutated genes increased in the resistant SCDCL (Fig. 3B). Mutations in ARID3B and PDS5A, both part of PyClone cluster 8, showed an inverse pattern with elevated cellular fraction in the sensitive SCDCL (2F4) and were absent from the resistant SCDCL (1G10; Fig. 3B). Mutations in typical colorectal cancer driver genes such as KRAS, TP53, and CTNNB1 were clonal and not associated with the respective response level of the SCDCLs.

In contrast to gene mutations, we did not detect any structural variants from the WES analysis that were enriched in either SCDCL. Copy-number variants (CNV) determined from WES data based on b-allele frequencies showed very similar gain and loss patterns in the tumor biopsy and the derived cell lines (Supplementary Fig. S3). However, the rather low purity of the tumor biopsy (≤40% purity) makes high-resolution CNV calls challenging. We identified 67 genes that differed in the absolute copy number between the resistant and sensitive SCDCL, and also showed significant differential expression (FDR < 0.1) in the same fold change direction as the copy number difference (Supplementary Data S7).

To decipher patterns of chromosomal heterogeneity on a single-cell level, we performed single-cell genomic analyses using multiplex interphase FISH (miFISH) of 14 gene loci within each nucleus (14). We designed a probe panel targeting the colorectal cancer genes COX2 (1q31.1), PIK3CA (3q26.32), APC (5q22.2), CLIC1 (6p21.32), EGFR (7p11.2), MYC (8q24.21), CCND1 (11q13.3), CDX2 (13q12.2), CDH1 (16q22.1), HER2 (7q12), TP53 (17p13.1), SMAD7 (18q21.1), SMAD4 (18q21.2), and ZNF217 (20q13.2).

The aneuploid cell population of the tumor biopsy consisted of a major clone in 39.8% of cells, and three minor clones (7.9%, 2.8%, and 2.8%) that showed copy-number alterations (CNA) of one gene from the major clone. The remaining cells (46.9%) were very heterogeneous. The chromosomal instability index (CII) was 43.3. The gain and loss patterns and the trajectory of the evolution of cell populations are shown in Fig. 4A.

The analysis of the derived parental cell population revealed a major clone in 54.6% of the population, which differed in one gene from the major clone of the biopsy. PIK3CA had no CNA in comparison with the major clone detected in the tumor biopsy, which had a loss of PIK3CA. We additionally detected seven minor clones ranging from 2.4% to 5.6% of the population that differed in one gene from the major population, suggesting a single chromosome segregation error. The remaining cells (20.2%) were dominated by a nonclonal distribution of genomic imbalances (Fig. 4B). The CII was 19.8. The rectal cancer specific distribution of genomic imbalances of the tumor biopsy was maintained in the derived cell population and in the SCDCLs. In addition, we applied a dedicated algorithm developed for the reconstruction of phylogenetic trees from miFISH data (33). The FISHtrees reveal that the resistant SCDCL 1G10, the sensitive SCDCL 2F4, as well as the parental cell population have a related evolutionary trajectory, yet develop clones specific to each of the cell lines (Supplementary Fig. S4).

Single-cell genetic analysis of the most resistant SCDCL (1G10) and the most sensitive SCDCL (2F4; Fig. 4C and D) revealed that the major clone of the parental cell population (54.6%) was present in an equal percentage (55.2%) in the resistant SCDCL, but in 81% of the cells of the sensitive SCDCL. Of note, the major clone of the tumor biopsy (39.8%) was observed in 4.6% of the cells in the parental cell population, expanded to 24.8% of the cells in the resistant SCDCL, but was absent in the sensitive SCDCL. This clone differed from the major clone by a loss of one copy of PIK3CA. This change could not be detected by WES because an approximately 20% shift in cellularity for a single copy change is beyond the resolution of WES (Supplementary Fig. S5). The CII was lower in the sensitive SCDCL (12.58) compared with the resistant SCDCL (19.39). Our results are consistent with the interpretation that the evolution of SCDCLs, despite ongoing chromosomal instability, is dominated by the selection of rectal cancer–specific genomic imbalances.

We subsequently analyzed gene expression profiles in the parental cell population, the three most resistant SCDCLs (1G10, 2D2, and 2G4), and in the three most sensitive SCDCLs (2F4, 1D2, and 2F1) by tRNA-seq (Supplementary Data S8). Global gene expression profiles separated the parental cell population and the resistant and the sensitive SCDCLs, displayed as a principal component analysis (PCA) in Fig. 5A. The differences in gene expression profiles between the three resistant SCDCLs and the pooled three most sensitive SCDCLs compared with the parental population are visualized as Volcano plots in Supplementary Fig. S6. To determine the functional space of these differentially expressed genes, we conducted pathway enrichment analyses using Ingenuity Pathway Analysis (IPA) and ToppFun (ToppGene tool suite; ref. 19). We found a significant enrichment for genes involved in Wnt-mediated β-catenin signaling, PI3K/AKT signaling, and Hedgehog signaling in the resistant SCDCLs (Fig. 5B). Activation z-score analysis implemented in IPA was used to identify pathways with differential activity between the resistant and the sensitive SCDCLs, relative to the parental cell population (e.g., candidate pathways will have increased activation in the resistant SCDCLs relative to the parental cell population, but reduced activation in the sensitive SCDCLs relative to the parental, or vice versa). This analysis confirmed both Wnt-mediated β-catenin signaling and PI3K/AKT signaling as having increased inferred activity in all three resistant SCDCLs relative to the parental cell population, but reduced inferred activity in the sensitive SCDCLs relative to the parental cell population. TCF7L2 was expressed highest in the most resistant SCDCL (1G10), and mRNA levels decrease with increasing sensitivity to radiation (Fig. 5C). The activation of the Wnt-mediated β-catenin signaling was confirmed independently using a Wnt reporter assay. As shown in Fig. 5D, SCDCL 1G10 has increased activity when compared with both the sensitive SCDCL 2F4 and the parental population (P = 0.05).

As an alternative to standard pathway enrichment, we also conducted preranked GSEA analysis using the MsigDB v6.2 canonical pathway collection. Genes were ranked using a t-statistic to produce a list sorted from the most significant upregulated gene down to the most significant downregulated gene for each paired comparison (each resistant SCDCL and the sensitive SCDCLs vs. the parental cell population). The Enrichment Map (21) plugin in Cytoscape v3.6.1 (22) was used to collapse pathways based on shared gene membership and functional overlap. For an initial comparison, we sought to identify any canonical pathway significantly upregulated or downregulated relative to the parental cell population. Clustered enrichment maps for the resistant SCDCL and sensitive SCDCLs are shown in Supplementary Fig. S7 illustrating the considerable differences in the significantly enriched pathways for each of these comparisons (red nodes are upregulated relative to the parental cell population, blue nodes are downregulated). For the most resistant SCDCL (1G10; Supplementary Fig. S7A), both Wnt-mediated β-catenin signaling and Hedgehog signaling were upregulated, while in the sensitive SCDCLs, there was clear downregulation of cell-cycle–related pathways, metabolic pathways, and DNA repair (Supplementary Fig. S7B).

Seven pathways exhibited normalized enrichment scores that were significant in at least one GSEA comparison, and correlated with the resistance phenotype (Fig. 5E). The most strongly correlated pathway is the Integrin-linked signaling pathway from the Pathway Interaction Database. This pathway plays a prominent role in cancer progression, and is implicated in all of the pathways described above (Wnt, Hedgehog, and PI3K/AKT signaling) that are more active in the resistant versus the sensitive SCDCL. Integrin signaling is also involved in mitotic spindle organization (34), similarly to ZWINT, which is mutated in the resistant SCDCL.

In summary, we established SCDCLs from a treatment-naïve rectal cancer biopsy. We found profound differences between individual SCDCLs in response to CRT. Genome and transcriptome profiling identified distinct gene mutations and signaling pathways associated with response to treatment.

Recent studies using multiregion sequencing of colorectal cancers revealed a considerable degree of ITH (12, 35). These studies also suggest that ITH results in the activation of tumor-promoting signaling pathways in different regions of the tumors; arguably, that would require tailoring therapeutic interventions taking the heterogeneity of individual tumors into consideration. In fact, studies in breast cancer provided compelling evidence that the selective pressure imposed by specific therapies can result in the expansion of resistant subclones, or subclones that utilize different pathways for survival (36). These studies underline the importance of ITH in therapy response. However, the functional role of individual clones, originating from single cells, in therapeutic resistance cannot be elucidated by such approaches. Clevers and colleagues recently reported the characterization of colorectal cancer organoids derived from single cells, which showed substantial differences in IC₅₀ values for chemotherapeutic agents (37). We hypothesized that the establishment of SCDCLs from a primary rectal tumor would provide suitable in vitro models to decipher the role of individual clones vis a vis response to CRT. Analyzing individual SCDCLs could allow identifying gene mutations, chromosomal aberrations, and genomic imbalances as well as different gene expression profiles in resistant compared with sensitive clones that are responsible for different therapeutic outcomes. Targeting specific pathways that mediate resistance in distinct subclones of a tumor could increase the proportion of patients who respond favorably to preoperative CRT; this would improve overall prognosis.

We established SCDCLs from a treatment-naïve biopsy of a patient with locally advanced rectal cancer, reflecting the genetic composition of the unperturbed tumor. Efficient in vitro expansion of human primary rectal tissue remains technically challenging; therefore, murine xenografting is commonly used to propagate primary colorectal tumor tissue ex vivo (38). PDX mouse models show close resemblance to the aberration profile of primary tumors (39) and maintain their intratumoral clonal heterogeneity (40, 41). We acknowledge that our approach provides only a snapshot of the ITH in the primary tumor. We did not attempt, nor would it be feasible, to reconstruct the complete clonal architecture of the primary tumor by SCDCLs. We realize that by establishing cell lines through PDX and subsequent single-cell cloning, a selective pressure is applied. The approach of establishing cell lines from a tumor biopsy, however, provides a model that allows to analyze heterogeneity from a functional point of view.

The model of SCDCLs allowed us to test the hypothesis that the genetic make-up of individual clones determines the degree of response to CRT of such distinct clones within a rectal carcinoma. Exposure of the SCDCLs to CRT in vitro revealed a gradient of different responses of individual SCDCLs. Of note, the SF of the parental cell population was intermediate, which likely indicates that the parental population is comprised of a mixture of clones with varying degrees of sensitivity to CRT. We conclude that our approach of establishing SCDCLs from a patient sample generates genuine in vitro models of response, and thus allows integration of function, phenotype, and genotype that could not be addressed by (multiregion) sequencing of tumor biopsies alone.

Using WES and miFISH, we could show that the clonal composition of the tumor biopsy was reflected in the derived parental cell population. Notably, we found that individual subclones present in the parental cell population were propagated in the SCDCLs, which is consistent with our previous results on single-cell cloning of colorectal cancer cell lines (42). We could recently show that the clonal composition of the parental population was in most cases reestablished after single-cell cloning maintaining specific genomic imbalances despite ongoing chromosomal instability. The clonal dynamics observed in our SCDCLs in this study is reminiscent of this phenomenon, suggesting that karyotype evolution is driven by the necessity to arrive at and maintain a specific plateau of chromosomal copy numbers as the drivers of carcinogenesis.

To further decipher genetic changes associated with treatment resistance, we performed WES-based clonal reconstruction. This revealed distinct patterns for the resistant versus sensitive SCDCLs with respect to subclones. Most striking was the enrichment of mutations in ZWINT and SV2A in the resistant SCDCL (cellular fraction of 49% and 46%) compared with both parental cell population (26% and 24%), tumor biopsy (both 2%) and sensitive SCDCL (both 0.9%), respectively. ZWINT is a kinetochore-associated protein that interacts with the kinetochore checkpoint protein ZW10 and is required for recruitment of ZW10 to the kinetochore (43). Silencing of ZWINT results in aberrant chromosome segregation (44). In addition, ZWINT expression has been implicated in resistance to multiple anticancer drugs (45). The potential functional role of missense mutation p.K74E remains elusive. However, the coincidence of this mutation and a higher CII in the parental cell population and resistant SCDCL, along with its enrichment in the resistant SCDCL might indicate functional relevance. Furthermore, GSEA identified multiple pathways involved in mitotic spindle formation as being upregulated in the resistant SCDCL and downregulated in the sensitive SCDCL, both relative to the parental cell population. Integrin-linked kinase signaling was the pathway most strongly correlated with the resistance phenotype; integrin signaling is required for accurate chromosome segregation (46). Synaptic vesicle protein 2A (SV2A) is a transmembrane protein that is critical for regulated exocytosis and neurotransmission (47). It is expressed in, for example, gastrointestinal stroma tumors (48). The potential role of SV2A with respect to sensitivity to CRT requires further exploration.

Mainly three signaling pathways were differentially regulated between the parental cell population and the resistant and sensitive SCDCLs: (i) Wnt, (ii) PI3K/AKT and (iii) Hedgehog signaling. Aberrant Wnt signaling is a defining feature of colorectal cancers (49). The Wnt signaling pathway was activated in the resistant SCDCLs as exemplified by pathway enrichment analysis and a subsequent Wnt reporter assay. The Wnt/β-catenin signaling pathway is tightly associated with resistance to CRT (3, 49). Enhanced Wnt/β-catenin signaling in the resistant SCDCL cannot be explained by the detected mutation in CTNNB1, which is clonal and not enriched in either the resistant or the sensitive SCDCL. However, we found that TCF7L2, one of the key transcription factors of the Wnt/β-catenin pathway, was overexpressed in resistant SCDCLs (50). We previously found overexpression of this transcription factor in patients that did not respond to CRT (7) and colorectal cancer cells were sensitized to CRT upon silencing of TCF7L2 (51). In the tumor from our patient, activation of Wnt/β-catenin signaling is most likely independent of the β-catenin mutation status. Besides, a heterogenous intratumoral distribution of β-catenin accumulated in the nucleus, a surrogate for Wnt/β-catenin pathway activation, has been described in colorectal cancers (52, 53). The second upregulated pathway in the more resistant SCDCLs was the PI3K/AKT signaling pathway, which is involved in cell metabolism, growth, proliferation, and survival (54). Folkvord and colleagues reported that patients with rectal cancer with a poor response to CRT had significantly elevated PI3K/AKT signaling measured by kinase activity compared with well-responding patients (55). Activation of PI3K/AKT signaling has been described as a major contributor to mediation of radioresistance in various cancers such as, for example, prostate and lung cancer (56, 57). Hedgehog signaling was the third pathway that was found to be upregulated in the more resistant SCDCLs. Hedghog pathway activation has been described to be associated with poor outcome after (chemo)-radiation in cervical cancer (58).

In conclusion, we show that patient-derived SCDCLs are representative models for deciphering the functional role of individual clones of a heterogenous rectal tumor in response to CRT. This approach allows us to identify specific gene mutations, chromosomal aberrations, gene expression profiles, and signaling pathways in resistant compared with sensitive SCDCLs that likely determine different therapeutic outcomes of heterogenous primary tumors. Future emphasis should be directed to the reconstruction of the clonal architecture of rectal cancers, which is the basis for specific targeting of pathways that mediate resistance in distinct subclones. This could increase the proportion of patients who respond favorably to preoperative CRT, improving overall prognosis.

No potential conflicts of interest were disclosed.

Conception and design: R. Braun, M.C. Cam, T. Ried

Development of methodology: R. Braun, Darawalee Wangsa, N.E. McNeil, K. Heselmeyer-Haddad, Z. Zhang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): R. Braun, L. Anthuber, D. Hirsch, Darawalee Wangsa, N.E. McNeil, K. Heselmeyer-Haddad, I. Torres, Danny Wangsa, M.A. Brown, A. Tubbs, P.R. Brauer, Z. Zhang, D.W. Rosenberg

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): R. Braun, L. Anthuber, D. Hirsch, Darawalee Wangsa, J. Lack, K. Heselmeyer-Haddad, I. Torres, Danny Wangsa, M.A. Brown, A. Tubbs, N. Auslander, E.M. Gertz, P.R. Brauer, D.L. Sackett, E. Ruppin, T. Ried

Writing, review, and/or revision of the manuscript: R. Braun, L. Anthuber, D. Hirsch, Darawalee Wangsa, J. Lack, N.E. McNeil, E.M. Gertz, D.L. Sackett, D.W. Rosenberg, T. Ried

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R. Braun, D. Hirsch, J.K. Habermann

Study supervision: R. Braun, A. Nussenzweig, E. Ruppin, T. Ried

The authors are grateful to Xiaolin Wu and Bao Tran (Frederick National Laboratory for Cancer Research, NCI, Frederick, MD) for sequencing. The authors thank Dr. Dali Zong and Dr. Alison McBride for help with irradiation of cells. The authors are grateful to Dr. Reinhard Ebner for critical comments, and Buddy Chen for editorial assistance. This research was supported, in part, by the Intramural Research Program of the NIH, NCI (Bethesda, MD). R. Braun and D. Hirsch were supported by Mildred Scheel postdoctoral scholarships of the German Cancer Aid (Deutsche Krebshilfe). L. Anthuber received a stipend of the University of Lübeck (Lübeck, Germany). D.W. Rosenberg was funded by NCI grant 5R01CA159976. We thank Jadwiga Jerman (Waterbury Hospital, Waterbury, CT) for clinical data and sample collection and Casey Dagnall (Cancer Genomics Research Laboratory, NCI, Gaithersburg, MD) for STR-profiling.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.