Monitoring Spontaneous Quiescence and Asynchronous Proliferation-Quiescence Decisions in Prostate Cancer Cells

The proliferation-quiescence decision is a dynamic process that remains incompletely understood. Live-cell imaging with fluorescent cell cycle sensors now allows us to visualize the dynamics of cell cycle transitions and has revealed that proliferation-quiescence decisions can be highly heterogeneous, even among clonal cell lines in culture. Under normal culture conditions, cells often spontaneously enter non-cycling G0 states of varying duration and depth. This also occurs in cancer cells and G0 entry in tumors may underlie tumor dormancy and issues with cancer recurrence. Here we show that a cell cycle indicator previously shown to indicate G0 upon serum starvation, mVenus-p27K-, can also be used to monitor spontaneous quiescence in untransformed and cancer cell lines. We find that the duration of spontaneous quiescence in untransformed and cancer cells is heterogeneous and that a portion of this heterogeneity results from asynchronous proliferation-quiescence decisions in pairs of daughters after mitosis, where one daughter cell enters or remains in temporary quiescence while the other does not. We find that cancer dormancy signals influence both entry into quiescence and asynchronous proliferation-quiescence decisions after mitosis. Finally, we show that spontaneously quiescent prostate cancer cells exhibit altered expression of components of the Hippo pathway and are enriched for the stem cell markers CD133 and CD44. This suggests a hypothesis that dormancy signals could promote cancer recurrence by increasing the proportion of quiescent tumor cells poised for cell cycle re-entry with stem cell characteristics in cancer.

To monitor and quantitatively measure the dynamic transitions of cell cycle states-from cytokinesis to S phase entry, we developed the Automated Temporal Tracking of Quiescence (ATCQ) analysis platform that includes a computational framework for automated cell segmentation (identification of individual cells in an image), tracking, cell cycle state identification, and quantification from movies (Supp Fig.1). The cell segmentation and tracking allows us to record the fluorescent reporter intensity changes within individual cells in real-time imaging, without the aid of a constitutive nuclear marker. The single-cell fluorescence changes over time, in turn, are used to obtain cell cycle state identification (G0, G1, or early S phase) and quantification, which allows us to examine the kinetics of the proliferation-quiescence transition.
The single-cell trace of fluorescent reporters, mVenus-p27Kand mCherry-hCdt1, graphed by ATCQ is consistent with what we manually observed in movies (Supp Fig.1).
We used the temporal changes in the fluorescence intensity readings of both reporters to create a classification that automatically assigned cell cycle phases for single cells (Supp Fig.1D). To quantitatively assess the relative change of two fluorescent reporters, we converted the relative change of the intensity readings from two adjacent time frames, t and t+1, into radian values (Supp Fig.1C). That is, the vector from time t to time t+1 on a plot of mCherry vs. mVenus fluorescence intensities is calculated into radian values over the range from -180° to 180°. A radial histogram of all recorded fluorescent vectors therefore appropriately reflects three primary clusters of cellular behaviors, corresponding to G0 entry, G0 exit or G1 entry, and G1 exit or early S entry. For example, a radian of -90° indicates cells exit G1 and enter into early S phase as mCherry signal decreases, while a radian of -180° indicates cells exit G0 and enter G1 as mVenus signal decreases while mCherry remains high.
This allows us to quantitatively visualize three distinct cell cycle behaviors at the proliferation-quiescence transition from the analysis of over 400 single-cell traces. The radial histogram exhibits the frequency and variability of three cell cycle behaviors in a population of asynchronously growing cells under normal, unperturbed conditions. The length of the bars depicts the frequency of individual cells with two fluorescent reporters exhibiting a specific behavior assigned to each radian. The clustering of cells into three distinct states is marked by the local minima of the frequency in radial histogram. Therefore, a given cell at any internal time is associated with state 1: G0, state 2: G1, or state 3: S phase. To facilitate the visualization of this, cells in G0 state are labeled in blue, while cells that exit G0 and progress into G1 phase are labeled in orange. When cells enter S phase, they are labeled in grey (Supp Fig.1D-F).
To assess the accuracy of our approach, we generated a plot to show all the cell-level radian assignment values for temporally adjacent movie frames. In this plot, all the radian values defining cell state can be compared among all cells in any movie frame t and the next frame t+1 (Supp. Fig. 1F). Each dot represents a single cell with the radian vectors drawn from adjacent time frames. Therefore, cells remaining in the same state within a time interval lie on the diagonal of the plot, while cells that transit between two states scatter above or below the diagonal. We considered cells to be correctly assigned if they either remained in their current state in the next frame, or cycle forward in the next frame. Our approach had a very high correct assignment rate, with 97-98% of frame intervals assigned with the correct cell-cycle transition order. This validates the feasibility of using ATCQ to investigate cell cycle dynamics.
Furthermore, because ATCQ uses the relative changes in intensity of the two cell cycle reporters over time to assign a cell cycle state, correct assignments are not affected by cell to cell variation in fluorescence intensity.

Validating ATCQ
Single-cell, real-time imaging produces large amounts of data, which is advantageous when assessing heterogeneous cell dynamics in a population. In order to quantify the proliferationquiescence decision with the large amount of imaging data, we used Kaplan-Meier (K-M) curves to estimate the time cells spent in a certain state (Supp. Fig.1H). The Kaplan-Meier estimator is a common statistical analysis tool to estimate clinical study survival functions that allows for individuals to be censored when they leave the study population (Kaplan and Meier, 1958). For our application, we took an advantage of it to estimate the probability of individual cells remaining in each state as a function of time. In the K-M curve, each step-down means cells exit that certain state. To compare time-in-states between treatment conditions, the time each cell spends in each state is first recorded, as well as whether a transition to the next state was observed. If the completion of the state was not observed, that cell was censored. Censorship would occur for example, if a cell underwent apoptosis or if it moved out of the observation frame.
To test whether ATCQ is able to capture impacts of environmental stimuli on the proliferation-quiescence decision, we performed a serum starvation to promote G0 entry. As expected, cells under serum deprivation exhibit a significant difference in G0 (blue) but not S phase. The increase in prolonged G0 upon serum starvation, is evident from the slower decline of the curve compared to the full serum conditions. However, cells also tended to have a longer G1 phase under serum starvation compared to full serum conditions. This result was initially unexpected, but may be consistent with a recent study showing that cells in G1phase can delay DNA replication in response to serum starvation (Cappell et al., 2016). ATCQ thus reveals dynamics of the proliferation-quiescence decision consistent with single-cell studies of others, and will facilitate further studies of this transition under various environmental perturbations (Supp. Fig.1H).
Previously, we revealed a role for PP2A complex in restricting Cdk2 activity after mitosis to modulate the proliferation-quiescence decision in tissue development (Sun and Buttitta, 2015).
We found that compromising PP2A activity caused 10% cells to fail to enter a developmentally regulated G0 state at the proper time. To examine whether PP2A affects G0 in cell culture, we performed a short-term (1h) treatment with the pan-PP2A inhibitor Okadaic Acid (OA) prior to imaging. As expected, analysis of the survival distributions showed that OA treatment results in a modest but statistically significant decrease in the median time in G0 (160 minutes vs. 140 minutes, p < 0.001). By contrast, no significant difference was observed for the time spent in the other two states (G1 and early S phase) (Supp Fig.I). Importantly the modest decrease observed in the G0 state is consistent with that expected if approximately 10% of cells are affected as previously shown (Naetar et al., 2014;Sun and Buttitta, 2015). This suggests that PP2A plays a role in promoting entry and maintenance of the G0 state in a fraction of cells and that ATCQ system can be sensitive enough to quantify the effects of perturbations that affect the proliferation-quiescence decision in a small fraction of cells.

Detailed Methods for ATCQ:
Live cell imaging: Cells were cultured at low density (to avoid contact inhibition) on 12-well plates in phenol red-free DMEM/10%FBS or 1%DMEM. Cells were imaged using an EVOS FL cell imaging system with a 20X objective lens or an IncuCyte Zoom at 37 o C, 5% CO2. The imaging intervals were 20-30 minutes. For PP2A inhibition, cells were exposed to Okadaic Acid (OA) at 50nM or DMSO (control) for 30 minutes. Media containing OA or DMSO was then removed, and fresh complete medium was added. Cells were imaged one hour later.

Image processing & analysis:
Individual cell identification (segmentation) and tracking were performed using a custom MATLAB program that automates these processes as described below. The program interface provides interactive adjustment of the parameters used at each step in the segmentation process, allowing the user to tailor the segmentation process to accommodate the characteristics of each image under analysis. Then, cell tracking is done with the Simple Tracker tool available on the MATLAB File Exchange. This tool uses the Hungarian algorithm to link cells from frame to frame into connected tracks using user-entered parameters for maximum movement distance allowed between frames and maximum gap in frames allowed between occurrences of the same cell track. Cells were tracked for the entire time they were visible in the field of view. A subset of tracks were verified by hand, to be sure that one track represented the same cell for the duration of the track. The output of this software is an Excel file containing the identified cells, the frame number(s) in which they appear, the size, approximate radius, and cell center location in each frame, and the average and maximum intensity values in the cell on both the red and green channels in each frame, as well as the segmentation settings used to obtain these results.

ATCQ Statistical analysis
Fluorescence Intensity smoothing: Analysis was restricted to cells tracked over at least 3 video frames. The mCherry and mVenus fluorescence signals were smoothed using a nearest neighbor proportion localized polynomial approach. The localized polynomial nearest neighbor value parameter was chosen using an Akaike Information Criteria (Akaike, 1974) based optimization with leave one out cross-validation. Calculations were performed using the Locfit package (2) in the software program R 3.0.1(R Core Team, 2014).

Transformation to survival time data:
For each cell, the cell cycle state over each pair of time points was assigned as described as before. For each cell, the aggregated time spent in each state was calculated. The event class for the last time point in each state was assigned as either "correct cycle" or "censored". Cells that were not assigned to progress forward in the cell cycle were censored from the analysis.
Statistical calculations: Associations between experimental time and overall fluorescence intensity were calculated via the Pearson's correlation coefficient and associated test (Tukey, 1949). Survival curves were generated via the Kaplan-Meier method and compared using the log-rank test. Survival time statistics were calculated using the survival package in R 3.0.1 (https://github.com/therneau/survival). Repeated measures linear regression was performed using generalized least squares regression using the nlme package with individual cells as the withinfactor level (Cox, 1972;Pinheiro et al., 2017).

Appendix 2: Du-145 Xenograft RNAseq comparison
Tumor cells were isolated by FACS from the CB17 SCID mice injected in the left cardiac ventricle with luciferase and GFP labeled Du145 PCa cells expressing a control shRNA construct from our previous report (Cackowski et al., 2017). Briefly, the presence of tumor was determined by bioluminescence imaging. Tibiae, femorae, pelvic and lumbar vertebrae bones with or without detectable tumor (involved vs. uninvolved) were dissected and crushed with a mortal and pestle to release marrow cells. Mouse cells were depleted with a mouse cell depletion kit (Miltenyi #130-104-694). Cells were incubated with an APC/Cy7 anti-human HLA A,B,C antibody (Biolegend #311426) and a PerCP/Cy5.5 conjugated anti-mouse hematopoietic lineage cocktail antibodies (BD Biosciences #561317) and DAPI. Cancer cells were isolated as single / viable / lineage negative / HLA + cells and collected in Qiagen RLT buffer containing βmercaptoethanol. Total RNA was isolated using the Qiagen RNeasy micro kit (#74004). After RNA quality assessment with an Agilent 2100 Bioanalyzer, we selected paired samples from two mice for further analysis. Libraries were prepared using the Takara Clontech SMART-seq V4 Ultra Low kit for single cell analysis -which is also appropriate for small numbers of bulk isolated cells. We then acquired approximately 40 million 50 nucleotide paired end reads per sample on an Illumina Hi-Seq 4000. Tuxeudo Suite software (Illumina) was used for analysis; including TopHat (version 2.0.13) and Bowtie2 (version 2.2.1) for alignment to the UCSC hg19 reference genome, FastQC for a second round of quality control, and Cufflinks (version 2.2.1) for differential expression analysis. Only sequences that aligned to the human (rather than mouse) genome were used for analysis. The 117 genes with a q value of < 0.05 for differential expression (involved vs. uninvolved) in either positive or negative direction were used for pathway analysis using DAVID version 6.8 software and the KEGG database. This dataset is publicly available in GEO under accession number GSE186061. All experimental procedures were approved by the University of Michigan Animal Care and Use Committee, protocol #PRO00007427.

Du-145 Xenograft RNAseq Results:
To discover new genes and pathways potentially involved in prostate cancer dormancy and recurrence, we compared gene expression from prostate cancer cells isolated macroscopic experimental bone tumors (actively growing) vs. prostate cancer cells isolated from bones of the same animals with no grossly detectable tumors (dormant). Tumors were detected by the presence or absence of bioluminescence signal from luciferase labeled human Du145 cells in SCID mice (Supp. Fig. 4) Samples and FACS isolation methods were from our prior study (Cackowski et al., 2017). We used DAVID version 6.8 software and the KEGG database to perform pathway analysis of genes differentially expressed in either the positive or negative direction and saw significant enrichment in the following pathways; arrhythmogenic right ventricular cardiomyopathy, TGF-beta signaling pathway, amyotrophic lateral sclerosis, basal cell carcinoma, and hippo signaling pathway (Supp. Fig. 4B). We interpreted the enrichment of TGF-beta signaling components as proof of concept for our approach as this pathway is known to be involved in regulation of prostate cancer dormancy (Bragado et al., 2013).

Flow Cytometry analysis and FACS
For flow cytometry analysis of DNA content, cells were fixed with 4% paraformaldehyde, and then stained with the FxCycle Violet DNA Stain (Thermo Fisher, F10347). Resuspended cells were analyzed on an Attune flow cytometer using the violet laser for DNA measurements. For cell sorting and assays in Supp. Figs. 2 and 4, cells were cultured in DMEM supplemented with either 10% FBS or 1%FBS for 24-72 hr, then subpopulations were sorted according to the intensity of their fluorescent reporters, using a BD FACS Aria II system. Cells were sorted into SDS-PAGE loading buffer or TRIZOL for immediate protein extraction or RNA isolation. A minimum of ~10 5 cells were collected for each experiment.

RT-qPCR
For RT-qPCR, RNA was isolated from sorted cells per the TRIzol manual, resuspended in water, and then treated with TURBO™ DNase to remove contaminating DNA. Using 500ng of RNA per sample, cDNA was synthesized using oligo(dT)20 and the Superscript III First-Strand Synthesis System (Invitrogen 18080051). qPCR using 0.5µL of cDNA per reaction was then