Tumor circadian clock strength influences metastatic potential and predicts patient prognosis in luminal A breast cancer

Abstract

Studies in shift workers and model organisms link circadian disruption to breast cancer. However, molecular circadian rhythms in noncancerous and cancerous human breast tissues and their clinical relevance are largely unknown. We reconstructed rhythms informatically, integrating locally collected, time-stamped biopsies with public datasets. For noncancerous breast tissue, inflammatory, epithelial-mesenchymal transition (EMT), and estrogen responsiveness pathways show circadian modulation. Among tumors, clock correlation analysis demonstrates subtype-specific changes in circadian organization. Luminal A organoids and informatic ordering of luminal A samples exhibit continued, albeit dampened and reprogrammed rhythms. However, CYCLOPS magnitude, a measure of global rhythm strength, varied widely among luminal A samples. Cycling of EMT pathway genes was markedly increased in high-magnitude luminal A tumors. Surprisingly, patients with high-magnitude tumors had reduced 5-y survival. Correspondingly, 3D luminal A cultures show reduced invasion following molecular clock disruption. This study links subtype-specific circadian disruption in breast cancer to EMT, metastatic potential, and prognosis.

Keywords: breast cancer; circadian data ordering; circadian medicine; metastasis; prognosis.

Fig. 1.

CYCLOPS reconstructed rhythms in noncancerous human breast samples. CYCLOPS was used to estimate sample circadian phase in noncancerous human breast tissue, including samples newly collected in Manchester, United Kingdom (N = 26), from GTEx (N = 167), and TCGA (N = 106). (A) Inferred time of peak transcript expression phase (acrophase) for select core clock genes in noncancerous human breast (outer, blue). The mouse paralogue acrophases (averaged across mouse tissues) are shown in the inner circle (inner, orange). In the mouse, time 0 was defined by the peak time of Arntl (Bmal1) expression. Rhythmicity was assessed by modified cosinor regression. All transcripts shown in blue had a regression P < 0.05. The human ordering was aligned to match the mouse acrophases. (B) Transcript expression was fit using cosinor regression including batch collection site. Batch-adjusted transcript expression is plotted as a function of CYCLOPS-predicted sample phase. The best-fit sinusoid and its significance are shown for each transcript. Solid lines represent BHq < 0.05. (C) Optimal alignment of CYCLOPS-predicted sample phases for the subset of time-stamped samples collected in Manchester (${\mathit{\text{Corr}}}_{\mathit{\text{Fisher}}}=0.69, p_{\mathit{\text{Fisher}}}<0.01$ ). (D) Histogram of inferred sample collection hour for noncancerous GTEx (filled) and TCGA (outlined) data. Sample phases are aligned as in (A). (E) Histogram of transcript acrophases including all significantly cycling transcripts in noncancerous breast tissues (BHq < 0.05). (F) PSEA was applied to CYCLOPS-ordered noncancerous human breast data. MSigDB hallmark gene sets with BHq < 0.05 are shown. Gene sets graphed furthest from the center had the most significant phase coordination. Gene sets highlighted in yellow also demonstrated phase coordination in mouse mammary tissues. (G) EnrichR was used to identify MSigDB Hallmark gene sets overrepresented among cycling genes in human noncancerous breast (BHq < 0.05, relative amplitude >0.33). The significance of pathway overrepresentation (−log BHq) is plotted against the significance of pathway phase coordination identified by PSEA (−log BHq). (H) Transcripts were ranked by cosinor regression F statistic. GSEA was applied to the ranked list. Hallmark gene sets with BHq < 0.1 are shown.

Fig. 2.

Subtype-specific circadian rhythm dysfunction in breast cancers. (A) Evaluation of core circadian organization in human breast cancer subtypes from TCGA data. Heatmaps depicting the Spearman correlation between selected core clock genes are shown for noncancerous breast tissue (N = 111), luminal A (N = 532), luminal B (N = 203), and basal/TNBC (N = 181). The zstat value and P-value are computed using a Mantel test and a reference correlation matrix of clock and clock-associated genes from the mouse atlas data. Higher Zstat scores denote a stronger resemblance to the established reference for healthy tissues. (B) Left, Representative images of mammary organoids derived from noncancerous breast and paired tumor tissues from the same individual. (Scale bar: 50 µm.) Right, Immunostaining demonstrates morphological structures. Hoechst 33342 was used for the nucleus and F-Actin was used for the cytoskeleton. (Scale bar: 50 µm.) (C) Circadian rhythms in matched normal and tumor organoids were monitored by a LV200 system. Representative bioluminescence images of organoids transduced with BMAL1-Luc from luminal A (N = 4) and TNBC (N = 3) subtypes were taken at 6-h intervals. (Scale bar: 100 µm.) Both raw and detrended signals are shown (red traces, tumor organoids; blue traces, nontumor organoids).

Fig. 3.

CYCLOPS reconstructed circadian rhythms from luminal A samples. CYCLOPS was used to estimate sample circadian phase in luminal A samples (N = 18 from Manchester and N = 193 from TCGA). (A) Acrophases are plotted for select core clock genes in luminal A tumors (outer, blue) and mouse paralogues (inner, orange). All transcripts in blue had a regression P < 0.05. The human ordering was aligned to match the mouse acrophases. (B) Transcript expression in luminal A samples was fit by cosinor regression including batch collection site. Separate models were fit to luminal A (red) and noncancerous samples (blue). Batch-adjusted transcript expression is plotted as a function of the CYCLOPS-predicted sample phase. Solid lines represent BHq < 0.05, and dashed lines represent a nonsignificant fit. (C) Histogram of luminal A transcript acrophases cycling with BHq < 0.05. (D) Histogram of the log amplitude ratio comparing luminal A and noncancerous genes, including transcripts that cycled in either luminal A or noncancerous samples with BHq < 0.05. (E) The estimated phase of luminal A tumor samples is plotted against the estimated phase of the matched noncancerous samples. (F) PSEA was applied to CYCLOPS-ordered luminal A data. Gene sets furthest from the center had the most significant phase coordination. All gene sets other than those highlighted in red also demonstrated phase coordination in noncancerous breast samples. (G and H) MSigDB Hallmark gene sets enriched for increased (G) and decreased (H) cycling in luminal A samples are shown. Transcripts significantly cycling in either luminal A or noncancerous samples (BHq < 0.05) were ranked by the log fold change in amplitude and analyzed by GSEA.

Fig. 4.

CYCLOPS clock magnitude varies in luminal A samples, correlates with EMT, and predicts prognosis. (A) Cartoon depicting CYCLOPS sample magnitude. CYLOPS projects the eigengene expression of each sample to a plane where the circular data structure is apparent. The angular position of any sample in the circle reflects its circadian phase. The radial distance from the center is the sample’s CYCLOPS magnitude (dotted line). This distance is a weighted sum of the amplitudes of individually cycling seed genes. (B) Histogram of luminal A CYCLOPS sample magnitudes colored by tertiles (N = 70). The distribution shows a long tail. (C) A histogram of transcript amplitude ratios comparing samples from the top-third and bottom-third magnitude groups. Transcripts that showed differential cycling by the CYCLOPS magnitude group were identified (BHq < 0.05). For differentially cycling transcripts, the cycling amplitude in the highest tertile was compared to the amplitude in the lowest tertile. (D) Batch normalized expression data are shown for representative core clock transcripts that exhibit differential cycling between magnitude groups (BHq < 0.05). Luminal A samples in the top tertile of CYCLOPS magnitude are shown in red and the bottom two tertiles in blue. The best-fit sinusoid is superimposed for each group. (E) The CYCLOPS estimated phase of luminal A tumor samples is plotted against the phase of matched noncancerous samples in TCGA. Samples are colored by the CYCLOPS magnitude group, distinguishing the top tertile from the bottom two tertiles. (F) Five-year patient mortality grouped by the tumor CYCLOPS magnitude group. Patient outcome data were obtained from the TCGA database. For each tertile/magnitude group, the percentage of patients who died within 5-y of diagnosis is shown. The increased risk in the high-magnitude tumor group remained statistically significant when evaluated in a logistic regression model that included patient age and the presence of known metastasis at diagnosis (P < 0.05). (G) Pathway overrepresentation (EnrichR) and enrichment analysis (GSEA) for transcripts that showed differential cycling (BHq < 0.05) among luminal A magnitude groups. The −log(BHq) from the EnrichR analysis is plotted against the normalized enrichment score (NES) from GSEA. (H) Batch normalized expression data for select transcripts in the epithelial–mesenchymal transition pathway. Data for the top tertile of CYCLOPS magnitude (red) and bottom two tertiles (blue) were fit to sinusoids.

Fig. 5.

Invasion of MCF-7 cells into the 3D collagen matrix is suppressed by circadian clock disruption. (A) Representative traces of circadian rhythms in BMAL1-Luc expression in MCF-7 cells transduced with lentiviral shScramble or shBMAL1. Left, raw data; right, normalized data. N = 3. (B) Hanging drop cell invasion assay of MCF-7 cells in a 3D collagen matrix. Representative images of cell migration are displayed at day 0 (Left) and day 3 (Right). (Scale bar: 100 μm.) N = 4. (C) The relative distance of cell migration was quantified and plotted. Data were normalized to WT which was set as 1. Unpaired t test, *P < 0.05, N = 4. (D) Representative traces of circadian rhythms in BMAL1-Luc expression in MCF-7 cells treated with KL001. N = 3. (E) Representative images of cell invasion assay of MCF-7 cells treated with KL001. (Scale bar: 100 μm.) N = 3. (F) Quantification of cell migration distance in (E), unpaired t test, *P < 0.05, N = 3.