Deciphering a global source of non-genetic heterogeneity in cancer cells

Abstract Cell-to-cell variability within a clonal population, also known as non-genetic heterogeneity, has created significant challenges for intervening with diseases such as cancer. While non-genetic heterogeneity can arise from the variability in the expression of specific genes, it remains largely unclear whether and how clonal cells could be heterogeneous in the expression of the entire transcriptome. Here, we showed that gene transcriptional activity is globally modulated in individual cancer cells, leading to non-genetic heterogeneity in the global transcription rate. Such heterogeneity contributes to cell-to-cell variability in transcriptome size and displays both dynamic and static characteristics, with the global transcription rate temporally modulated in a cell-cycle-coupled manner and the time-averaged rate being distinct between cells and heritable across generations. Additional evidence indicated the role of ATP metabolism in this heterogeneity, and suggested its implication in intrinsic cancer drug tolerance. Collectively, our work shed light on the mode, mechanism, and implication of a global but often hidden source of non-genetic heterogeneity.


INTRODUCTION
It has been well established that, although cells in a clonal population share the same genetic background, they are often highl y hetero geneous, with genes exhibiting significantl y variab le e xpression le v els (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11).Heterogeneous gene e xpression can lead to phenotypic variations at the single-cell le v el and has been studied in di v erse systems and organisms, fr om micr obial str ess r esponse to cancer progr ession, wher e it plays critical biological functions and has important clinical relevance (12)(13)(14)(15)(16)(17).Thus, understanding the mechanism and function of non-genetic heterogeneity is critical for precisely intervening with diseased cells.
Studies on the source of non-genetic heterogeneity have been largely focused on individual genes ( 14 , 18-21 ), wher eby the expr ession variabilities of specific genes can confer phenotypic heterogeneities that have significant biological implications.For example, heterogeneous expression of a stress protectant allows bet-hedging in microbial cells ( 22 ), variability in target genes of pluripotency regulators may underlie de v elopmental priming in early mouse embryos ( 23 ), and variable expression level of the secretory leukocyte peptidase inhibitor gene affects cancer clonal dominance ( 24 ).Notably, in these examples, the preexisting, rather than induced, non-genetic heterogeneity can play unexpected roles in affecting the responses to perturbations, highlighting the biological significance of intrinsic clonal variability.
In addition to non-genetic heterogeneity arising from the expression variability of individual genes, recent advances in genomics have implicated the presence of cell-to-cell variability in the expression of the entire transcriptome, as specific perturbations such as transcription factor ov ere xpression or drug treatment can alter the size of the transcriptome by globally modulating gene transcription (25)(26)(27)(28)(29).Such global modulation is distinct from gene-specific modulation as it would globally upregulate or downregulate the expression of all genes nonspecifically in a cell.Global modulation of gene transcription has also been inferred by single-cell assays characterizing transcript le v els or transcriptional dynamics of multiple genes in the same cell (30)(31)(32)(33).Howe v er, because of the technical challenges for detecting and analyzing globally modulated transcription in time, its contribution to non-genetic heterogeneity often remains hidden ( 32 ).Mechanisticall y, w hile cell size ( 34 , 35 ) and mitochondrial le v el ( 36 , 37 ) hav e been implicated in the control of global transcription rate, how such factors temporally modulate global transcription in single cells remains largely unclear.Functionall y, w hile it has been implica ted tha t the modulation of transcriptome size can be induced to confer enhanced proliferation or chemotherapeutic resistance for cancer cells ( 25 , 26 , 29 ), the role of pre-existing cell-tocell variability in the expression of the entire transcriptome remains elusi v e.
Temporal analysis of noise or fluctuations leading to nongenetic heterogeneity is necessary for a systematic understanding, as the temporal characteristics of fluctuations can play a pivotal role in determining the mode and function of the resulting clonal heterogeneity ( 38 ).As an example, temporal analysis of individual genes' expression fluctuations has r ecently r e v ealed that the e xpression status of a small fraction of genes can be heritab le ov er generations in multiple cancer cell types ( 39 , 40 ), leading to long-lasting (and static) non-genetic heterogeneity in a clonal population and potentially conferring pivotal cancer phenotypes such as drug resistance and metastasis ( 39 , 40 ).Consequently, it would be highly desirable to carry out temporal analysis of fluctuations influencing the expression of the entire transcriptome.
In this work, we set out to temporally analyze the global modulation of gene transcription and decipher the resulting non-genetic heterogeneity in the expression of the entire transcriptome.To probe global transcriptional modulation, we devised an approach based on multi-integrated r eal-time transcriptional r eporters, allowing us to temporall y anal yze modulated global tr anscription r ate.We found that the instantaneous global tr anscription r ate (i.e. the ra te a t some instant in time) exhibits temporal fluctuations at the cell cycle timescale, whereas the time-averaged rate is variable among cells, indicating the co-existence of dynamic and static modes of heterogeneity.These observations were corroborated with analyses using public single-cell transcriptome data.Mechanistically, a cell-cycledependent ATP metabolic switch likely contributes to dynamic hetero geneity, w hile heritable differences in the timeaveraged metabolic state among cells could give rise to static hetero geneity.Functionall y, the heritable global transcription rate was implicated in differential intrinsic drug tolerance in cancer cells.Altogether, our temporal analysis of the global transcription rate yielded key insights into a global source of non-genetic heterogeneity, opening up avenues for broadl y anal yzing cell-to-cell variability at the global scale.

Plasmid construction and virus packaging
Molecular cloning was performed according to standard protocols, and plasmids were replicated in DH5 ␣ Esc heric hia coli .For the plasmid used for the multi-integrated reporter system, TRE3G promoter sequence, CFP sequence, and 24 × PP7 sequence were obtained from existing plasmids in the lab and then inserted into the piggy-Bac backbone plasmid.For plasmids used for two-color reporter system, endogenous promoters of ACTB ( −1284 bp to +1052 bp) and UBC ( −334 bp to +877 bp) were PCR amplified from U2OS genomic DNA and then assembled into 24 × MS2 plasmid by Gibson assembly.For cell cycle reporter NLS-Citrine-PCNA, the PCNA cassette was constructed by assembly PCR from synthetic DNA oligos, which was then Gibson assembled into lentivirus transfer backbone together with CMV promoter and Citrine cassettes.For A TP biosensor, the A TP sensor sequence was obtained from iATPSnFR1.0(Addgene #102550), and SV40 NLS sequence was added to the C-terminus to ensure the nuclear localization of the A TP sensor.A TP sensor cassette was synthesized (RUIBO) and inserted into the lentivirus transfer backbone by enzyme ligation.For plasmid validation, 24 × PP7 and 24 × MS2 repeat sequences were confirmed by enzyme digestion, and all the plasmid sequences were confirmed by Sanger sequencing (RUIBO).
Virus packaging was performed using the threeplasmid system.HEK293T cells at 60-80% plating density were transfected with lentivirus plasmids using polyethyleneimine (PEI).Medium was changed 8 h posttransfection and the culture supernatant containing viral particles was collected two days post-transfection.

Simulation-assisted reporter design
We simulated the transcriptional dynamics based on a two-state model using the Gillespie algorithm.Genes turn on / off randomly based on the thermodynamic distribution, and transcription only occurs during the on state: The parameters used in our simulations are: where k on , k of f , k m were estimated from our imaging experiment of 3-B5-6 cell line in the condition of 0.012 ng / L doxy cy cline, k s was tuned so that the total transcription le v el is comparab le to the e xperimental results.To model the effect of global fluctuations, we multiplied the parameter k m with a temporal signal r epr esenting the global fluctuations, for example, a sine wave function for periodic fluctuations or smoothed random time series for random fluctuations.The rationale for choosing to modulate k m is that RN A pol ymerase elongation can be globally regulated in the cell ( 37 ).The amplitude of global fluctuations ranges from 10% to 80%, and the period ranges from 4 to 48 h.The rationale for choosing such an amplitude range or period range was based on experimental measur ements (Figur e 3 H, Supplementary Figure S4D).For a single gene, its transcriptional activity dynamics are quantified as the number of nascent mRNAs being produced over time.Total simulation time lasts 48 h and the sampling interval is 10 min.
The minimum number of reporters (i.e. the number of reporters is designated as n) needed to quantify global fluctuations is measured as the minimal n when the CV 2 of the sum of n traces (i.e. each reporter gi v es rise to one nascent transcript count trace) is less than twice the CV 2 of global fluctua tions, so tha t the magnitude of the global fluctuations is comparable to the intrinsic noise.The rationale for such a criterion is multi-faceted.First, based on experimental results, intrinsic fluctuations are stronger than global fluctuations (Supplementary Figure S4B), i.e.CV i ntri nsi c > CV global .Thus, it is necessary to reduce intrinsic noise by av eraging ov er multiple traces until the strength of intrinsic noise is comparable to that of global fluctuations.Second, it should be noted that the choice of the criterion, i.e.
is empirical.This is because only when the intrinsic noise strength is reduced (by averaging) to a level that is comparable to global fluctuations, one can unambiguously detect the temporal waveform of global fluctuations.Please refer to Note S1 for more details.
In Supplementary Figure S3, we simula ted dif fer ent r eporter design schemes, and reporters from different regulons were assigned with different promoter switching rate constants, i.e. for promoter X, both k on and k of f are onefifth of the values above and for promoter Y, both parameters are fiv e-fold of the values above.

Stable cell line construction
To construct the multi-integrated reporter cell line, plasmids containing TRE3G-CFP-24xPP7 and rtTA were cotransfected into U2OS cells with the piggyBac transposon plasmid by Lipofectamine LTX (Thermo).Single cells containing stably integrated reporter plasmids were FACSsorted based on the CFP signal into individual culture wells.Many monoclones were infected by lentivirus containing PCP-mCherry and were then screened for the number of nascent reporter transcriptional sites by microscopy.A monoclone with an appropriate number of reporter transcriptional sites was identified and named 3-B5 (third plate, B5 well).To minimize the influence of the heterogenous expression of PCP-mCherry, a monoclone deri v ed from 3-B5 cells was chosen and named 3-B5-6, which was the reporter cell line used in this wor k.To enab le the quantification of cell cycle phase, the 3-B5-6 cell line was infected by lentivirus containing NLS-citrine-PCNA, and the resulting cells were FA CS sorted.To ena ble the quantification of ATP le v el, the 3-B5-6 cell line was infected by lentivirus containing the aforementioned ATP biosensor, and the resulting cells were FACS sorted.For the construction of dual PP7 / MS2 reporter cell line, U2OS cells co-expressing PCP-mCherry and MCP-CFP were FACS sorted after lentivirus infection, and the resulting cell line was transfected with four different plasmids containing TRE3G-YFP-24xPP7, rtTA, endogenous promoter-iRFP-24 × MS2 and piggyBac transposase by Lipofectamine LTX.Cells expressing four differ ent fluor escent proteins (mCherry, CFP, iRFP and YFP) were FACS sorted.

Identification of reporter integration sites
To identify the integration sites of the reporter cell line 3-B5-6, we performed whole genomic DNA sequencing.Genomic DNA of 3-B5-6 was extracted by using a universal genomic DNA kit (CWBio, CW2298).Whole genomic sequencing was performed by GENWIZ.Indels of human genome were mapped by pindel, which yields plasmid integration sites.

Time-lapse fluorescence microscopy
U2OS cells were seeded into a 24-well glass-bottom plate (Cellvis) 48 h prior to imaging.The culture medium was replaced with phenol r ed-fr ee medium 8-12 h prior to imaging.Time-lapse images were acquired with an automated inverted microscope Nikon Ti-E with a Plan-Apo 40 ×/ 0.85 objecti v e (Nikon).The microscope was equipped with a white-light LED (Lumencor SOLA) and an sCMOS camera (Hamamatsu ORCA-Flash4.0V2).Standard fluorescence microscopy filter cubes (Chroma and Semrock) were used for the acquisition of CFP , GFP , YFP , mCherry and iRFP fluorescence.The sample stage was maintained under 37 • C and 5% CO 2 humidified air in a custom environmental chamber.Multi-channel fluorescence images were acquired e v ery 10 min using a custom Micro-Manager ( 41 ) (v1.4) script for 24 to 48 h.
For imaging experiments with the cell line 3-B5-6 / Citrine-PCNA, doxy cy cline was added 12 h prior to imaging to ensure stead y-sta te transcription.The medium was replaced with fresh phenol red-free medium with doxy cy cline of the same concentration before imaging.For cell cycle blocking experiments, 2 mM thymidine (Sigma-Aldrich, T9250) was added after ima ging.Ima ges of mCherry, YPF, CFP and brightfield were captured every 10 min for 48 h.For imaging nascent transcription sites, four z-stacks images with 1.2 m spacing wer e acquir ed to capture loci in different focal planes.

Image analysis
Fluor escence images wer e analyzed using a custom MAT-LAB GUI program that allows automatic analysis as well as manual correction.Images were first background corrected, and segmented to obtain cell masks, which were then tracked across frames to obtain cell trajectories.More specifically, for background correction, the image was fitted using a 2D second-order pol ynomial function, w hich was used to remove the background and flatten the original image.The resulting image containing fluorescentlylabelled cell nuclei was edge detected using the Sobel opera tor, dila ted, and inner filled to obtain nuclei masks, which were then tracked based on the distances between adjacent fr ames.The tr acks of individual n uclei were man ually corrected with the GUI.
For nascent transcription site detection, the z-stacks images were firstly background eliminated, maximum intensity projected, and filtered by a Laplacian of Gaussian filter.The pixel positions of local maxima that pass a predefined threshold were regarded as nascent transcription sites.Traces of single loci along time were tracked based on the distance between adjacent frames and were corrected manually.The intensity of individual loci was calculated as sum(signal-bg) / bg, where the signal was defined as the in-tensity of a 3 × 3-pixel square around the loci pixel and the background (bg) was defined as the median intensity of a 9 × 9-pixel square with a 5 × 5-pixel hollow around the loci pixel.When the loci were below the detection threshold, their intensities would be set to zero.Global transcription ra tes a t specific time points were calcula ted as the summed intensities of all detected transcription sites in a cell.Mean transcriptional intensities or the number of activated transcription sites were also extracted.Autocorrelation analysis was performed on each cell's global tr anscription r ate tr ajectory, and the time lag of the first maximum was used as an estimate of the period of the dynamics.

Cell cycle phase measurement
The cell cycle phase was determined using both the Citrine-PCNA signal and morphological characteristics of the nucleus.More specifically, M phase was identified by the disappearance of nuclear boundaries.During S phase, nuclear Citrine-PCN A aggregates, w hich can be characterized as the emergence of bright fluorescence puncta or spatially une v en nuclear fluorescence intensities.During the S to G2 transition, Citrine-PCNA gathers to specific genome positions and then dissipates quickly.Because of these characteristics, we defined a parameter PCNA var based on the   variance of the PCNA signal as follows: PCN A var = max ( I YF P ) − mean ( I YF P ) mean ( I YF P ) .
Thus, G1 / S transition is defined as the first rise of PCNA var , and S / G2 transition is defined as the position of max value of PCNA var after G1 / S transition.

Calculation of intrinsic and extrinsic noise
Gene-specific and gene-nonspecific fluctuations (analogous to intrinsic and extrinsic noise) inside a cell were calculated based on the method previously introduced ( 2 , 42 ).Note that we aimed to compare the strengths of fluctuations within a single cell using temporal dynamics of multiple gene loci, which is distinct from typical operations using static snapshots of a population of cells.Because there are more than two dynamic traces in a cell in our experiments, we enumerated all pair combinations in a single cell.For example, if a cell has 10 traces, then 10 × 9 / 2 = 45 pairs would be enumerated.The intrinsic (gene-specific) noise ( η int ), extrinsic (gene-nonspecific) noise ( η ext ), and total noise ( η tot ) of each pair is defined as follows: Here, E 1 and E 2 are the two tr anscription tr aces in the chosen pair.The three noise terms were then calculated for all pairs inside each cell and the root mean squares of all pairs were computed as a cell's fluctuation strengths.

Quantification of dynamic and static heterogeneity using temporal data
Temporal dynamics from our reporter system offer unique opportunities for systematically decomposing variability in the observed signal (i.e.global transcription rate) into static and dynamic components (related to Figure 3 H).More specifically, if we define the global transcription rate at a specific time point (t) in a specific cell (c) as E t,c (i.e.instantaneous global tr anscription r ate), then the total heterogeneity of the population can be defined as: Mean ( E t,c ) by t,c .
To quantify dynamic heterogeneity, we first average rate le v el across cells and then compute the variability along the cell cycle: ) by c by t

Mean Mean ( E t,c ) by c by t
.
Operationally, to compute the variability along the cell cycle, we need to interpolate the global transcription rate trajectory according to the cell cycle phase measured by the PCNA reporter (i.e. between 0 to 100 percent) in order to average across cells.
To quantify static heterogeneity, we first compute the time-aver aged r ate le v el in indi vidual cells, i.e.E c = Mean ( E t,c ) by t .Static heterogeneity can then be calculated as: Mean ( E c ) by c .
The residual heterogeneity is thus: It should be noted that there are other ways to decompose the heterogeneity, for example, CV dynamic and CV residual could be merged into a single term when CV 2  dynamic is defined as the weighted average of individual cells' CV 2  dynamic .Yet, we adapted the above decomposition method because we focused on extracting and comparing CV st at ic and CV dynamic , without being influenced by the intrinsic noise in the summed transcriptional activity signal.In other words, in the current decomposition method, the intrinsic noise of the 'pseudo reporter' (i.e. the sum of all reporter alleles) is accounted for by the last term, CV residual .To determine if static hetero geneity simpl y arises from random fluctuations in gene transcription, we assume that the number of acti v e transcription sites observed during a defined time window within a cell obeys a binomial distribution E c random ∼ B( N, p ) .Here, N = 10 total gene loci in our system, and p can be calculated by the mean number of acti v e sites ov er total gene loci N .The binomial distribution is then multiplied by the mean transcriptional activity per gene loci to estimate the static heterogeneity arising from random fluctuations.

PEG-mediated fusion experiment
An entire well of 6-well U2OS cells was first concentrated into 25 l by centrifugation.The centrifuge tube was kept in a 37 • C bath, and 25 l 50% polyethylene glycol (PEG, MW 8000) was added drop by drop.The tube was then gently shaken to bring the PEG into contact with the cell and remained still for 90 s. 500 l pre-warmed 37 • C Opti-MEM (Gibco) was then added, and the sample tube was further kept in the 37 • C bath for 30 min.Cells were then centrifuged and washed with FBS-free culture medium twice and then seeded onto glass-bottom wells for imaging.After seeding, the medium was replaced with a phenol r ed-fr ee medium containing 0.01 ng / l doxy cy cline.Cells were imaged for 24 h with a frame rate of 10 min per frame.

ATP level quantification
To test whether ATP le v el correlates with nascent transcriptome or whether ATP le v el is heritable across generations, nuclear ATP le v els in indi vidual cells wer e measur ed using the aforementioned ATP biosensor.Nuclear ATP le v el was quantified as the mean value of GFP intensity inside the nucleus mask.

Perturbation of ATP metabolism
To perturb ATP metabolism, we used 2-Deoxy-D-glucose (2DG, Sigma-Aldrich, D8375) or oligomycin A (Selleck, S1478) to inhibit different stages of ATP metabolism.For experiments measuring transient responses to inhibitors, cells were prepared as in the time-lapse imaging section, and 5 M oligomycin or 10 mM 2DG was added after 5-10 h of ima ging, and ima ging continued for 24 h.For experiments measuring stead y-sta te responses to inhibitors, inhibitors were added 24 h prior to ima ging, and ima ging continued for 48-72 h.For experiments to measure the transient responses at the S cell cycle phase to inhibitors, 2 mM thymidine was added 24 h prior to imaging.

Nascent RNA labeling
Nascent RNA labeling was performed by EU (5-Ethynyl Uridine, Abcam ab146642) incorporation.Cells were incubated in the culture medium supplemented with 1 mM EU for 1 h, and were washed by DPBS before fixation in 4% paraformaldehyde.Cellular fluorescence intensity of the ATP sensor and positions of cells were recorded by microscop y befor e cool ethanol permeabilization.Then incorporated EU was labeled by Click-iT chemistry using Azide Flour 454 dye (Sigma, 760757) or Alexa 488 dye according to manufacturer`s instructions (Invitrogen, C10329).The fluorescence of labeled EU from the recoded sample positions was then measured on the microscope.

Imaging-based intrinsic drug tolerance experiments
To quantify the relationship between global transcription rate and intrinsic drug tolerance, we imaged the 3-B5-6 / Citrine-PCNA cell line under 0.01 ng / l doxy cy cline (added 12 h prior to imaging) for 24 h before adding indicated chemotherapeutic drug.After drug addition, imaging continued for 72-96 h at the rate of one frame per hour.The death of cells was identified manually by the disappearance of nuclear fluorescence signals as well as the detachment of cells from the surface.In the analysis, the average global tr anscription r ate over the first 24 h without drug was used to quantify the intrinsic global transcription rate in individual cells.

FACS-based intrinsic drug tolerance experiments
To deri v e an assa y f or analyzing intrinsic drug tolerance in a relati v ely high-throughput format, we need to have enough numbers of cells with either low or high global transcription rates and then subject them to drug treatments, followed by cell survival analysis.We thus resorted to flow sorting of CFP that is expressed under the control of doxycycline induction in the 3-B5-6 / Citrine-PCNA cell line.We quantified CFP induction (0.01 ng / L doxy cy cline) ov er a 1-day time course and correlated the intensity of CFP with integrated global transcription rates.By doing so, we found that the two quantities exhibited a maximal correla tion a t about 12 h post-induction.We concluded that CFP intensity at this time point could be used as a proxy for global transcription rate, and thus sorted based on CFP to obtain cells with either low or high global transcription rates (corresponding to the lowest or highest 20% of CFP intensity).Sorted cells were then plated on 96-well glass bottom plates, allowed for ∼12 h for cell attachment, and drugs were then added at indica ted concentra tions.Following drug addition, images were taken for at least 12 fields of view to analyze cell numbers, and survival cell fraction at 24, 60 or 72-h post drug addition (depending on the lethality of the drug) was used to compare drug tolerance of cells with low or high CFP.
We took 29 drugs from the Tocriscreen Epigenetics Library (Tocris , Cat.No .6801) that hav e been pre viously used in the literature as tumor-killing drugs (Supplementary Table S1).Out of these drugs, 11 drugs did not display cell killing within three days, one drug caused immediate lethal effect (pre v enting useful quantification of tolerance), and the remaining 17 drugs displayed cell killing effect within three days (Supplementary Table S1).

Bioinformatic analysis
For the analysis of single-cell RNA-seq data of U2OS cells, data from GSE146773 and https: //github.com/CellProfiling/SingleCellPr oteogenomics were downloaded, which contained intronic and exonic read counts as well as spiked-in ERCC read counts.We retained cells with ERCC read counts between 100 and 600, and removed genes detected in less than 100 cells.To quantify how the global transcription rate is modulated across the cell cycle, U2OS cells were divided into 14 bins based on the pr e-sequenced FUCCI r eporter signals, and summed spiked-in normalized total intron read counts were calculated for cells in each cell cycle bin for three different sets of genes, i.e. low expression genes (bottom 30%), medium expression genes (30%-70%), and high expression genes (top 30%).To infer potential regulators of global transcription rate, we used summed spiked-in normalized total intron read counts in individual cells as a 'target gene', and used all genes' normalized exon read counts as 'regulators' for inputs of GEINE3 (v1.14.0).Gene ontolo gy anal ysis of the top 100 genes from GEINE3 was performed by clusterProfiler (v4.0.2).For the analysis of single-cell RNA-seq data of A2780 cells, raw sequencing data from Wang, 2021 (GSE162256) was downloaded and processed.To map the A2780 data, we used Trim Galore (v0.5.0) to remove adaptors and used STAR (v2.5.3a) to map to STAR index generated using spiked-in sequence containing hg38 GENCODE v34 genome.TPM calculator (v0.0.3) was used to compute intronic and exonic read counts, which were then normalized by individual cells' spiked-in read counts.Because three spiked-in sequences in this dataset varied by 100 times in read counts, we used the spiked-in sequence with the highest expression for normalization.A2780 cells with spiked-in reads larger than half a million or with lower than 10 5 read counts were filtered out, so were genes detected in less than 5 cells.
For the analysis of drug tolerance using cancer cell line datasets, both gene expression profiles (RNA-seq, 21Q3 release) and drug sensitivity data (Drug sensitivity AUC (Sanger GDSC2) and Drug sensitivity IC50 (Sanger GDSC2)) were downloaded from the Cancer Dependency Map w e bsite ( https://depmap.org/portal/do wnload/ ).A dataset containing 597 cancer cell lines and 153 drugs was obtained.Note that not all cell lines have drug sensitivity data for all drugs.Using the expression matrix for the 597 cell lines, we performed hierarchical clustering analysis of the top 100 genes from GEINE3 based on expression correlation, and used genes in the largest co-expression module ( n = 38 genes) to calculate the modulation index in each cell line, defined as the median expression level of this module.For individual drugs, correlation analysis between drug tolerance measure (AUC or IC50) and modulation index across cell lines was performed.

Analysis of public scRNA-seq data implicated the presence of non-genetic heterogeneity in the global transcription rate
While it has been implicated by recent single-cell transcriptome studies that the total amount of mRNAs in a single mammalian cell can be variable within a clonal cell population ( 43 , 44 ), it is unclear whether such non-genetic heterogeneity in transcriptome size r epr esents heterogeneity in the rate of global transcription or heterogeneity in the rate of global mRNA degradation, as the total mRNA abundance is determined by the rate of global transcription together with the rate of global mRNA degradation.Addressing this question would be critical for understanding the cause and consequence of non-genetic heterogeneity in transcriptome size (i.e. total mRNA abundance).
To explore this, we resorted to public single-cell transcriptome data of the human osteosarcoma U2OS cell line with spiked-in RNAs ( 45 ), using which we could better estimate the transcriptome size per cell compared to data without spiked-ins ( 43 , 44 , 46 ).We found that, indeed the estimated transcriptome size displays a large cell-to-cell variability (Figure 1 A; Materials and Methods).To investigate the potential cause of such variability, we first analyzed the rate of global transcription in individual cells by leveraging the intronic reads in RNA-seq data, which r epr esent unspliced pr e-matur e RNAs ( 47 , 48 ) (Supplementary Figure S1A).Importantly, intronic reads can be used as a proxy of the nascent RNA synthesis rate and have been widely used in the analysis of gene transcription kinetics (47)(48)(49)(50)(51).We thus used intronic read counts (spiked-in normalized) summing over all genes in each cell as a proxy of the rate of global transcription, and found that the total intronic reads yield a high explanatory power for the transcriptome size (i.e. total exon read count) (Figure 1 B, left).We next estimated the rate of global mRNA degradation in individual cells using the ratio of total intronic read count over total exonic read count in a cell (Supplementary Figure S1A), analogous to previous gene-level estimations ( 47 , 51 ).It should be noted that such estimations are built on steady-state assumptions (Supplementary Figure S1A).We found that this rate yields a much weaker explanatory power compared to the global transcription rate (Figure 1 B, right).Similar observations were made for another cell line (Supplementary Figure S1B-D).These results indicate that compared to heterogeneity in the global degradation rate, the heterogeneity in the global transcription rate can better explain the heterogeneity in transcriptome size.
Thus, it appears that the global tr anscription r ate is highl y hetero geneous across cells, w hich is an often-hidden source of non-genetic heterogeneity for mammalian cells, as most studies focused on cell-to-cell variability at the le v el of specific genes.Many aspects regarding the heterogeneity in global tr anscription r ate remain elusi v e; for e xample, w hether such hetero geneity is dynamic (i.e.hetero geneity arising from dynamic fluctuations) or static (i.e.heterogeneity arising from co-existence of cells with different states), how it is modula ted, and wha t relevance it might have (Figure 1 C, Supplementary Figure S1E).

Simulations guided the design of a reporter system that enables dynamic analysis of global transcription rate
To address these questions, it would be necessary to temporall y anal yze the global tr anscription r ate in individual cells, because snapshot-based methods such as singlecell sequencing cannot distinguish between dynamic and static modes of heterogeneity.We thus aimed to de v elop an imaging-based reporter system to analyze global transcription rate dynamics in single cells.
Analyzing global transcription ra te d ynamics r equir es capturing instantaneous global transcription rates.Howe v er, it is almost impossible to quantify time-lapse global tr anscription r ate by measuring all genes in a single cell.We thus proposed to measure a proxy of global transcription rate by using a reporter system based on multi-integrated transcriptional r eporters (Figur e 2 A top) (52)(53)(54).These r eporters can report real-time transcriptional activity because the RNA stem-loops on the nascent transcripts would recruit GFP or mCherry-labeled RNA binding proteins such that nascent transcription sites would appear as bright fluorescent foci due to the local enrichment of RNA binding proteins.
To illustrate the rationale behind such an approach, we needed to first discuss the stochasticity involved in gene tr anscription.The tr anscriptional dynamics of a gene are typically described by a two-state model whereby the promoter switches stochastically between ON and OFF states, and RNAs are synthesized during the ON state, leading to stochastic bursts of RNA production (Supplementary Figure S2A).The stochasticity in transcriptional dynamics can arise from gene-specific fluctuations (such as the noise inherent to the transcriptional reactions), as well as from gene-nonspecific fluctuations (such as fluctuations in the cellular sta te tha t af fect many genes) ( 2 , 42 ).Because gene-specific fluctuations are temporally uncorrelated between genes while gene-nonspecific fluctuations ar e shar ed across genes, summing transcriptional activity dynamics from multiple genes could in theory, remove gene-specific fluctuations while maintaining gene-nonspecific fluctuations.And as illustrated in the cartoon, the multi-genesummed activity can serve as a proxy of the global transcription rate, allowing capturing the global modulation of transcriptional activities (Supplementary Figure S2A).
We next used computer simulations to validate the multireporter approach and to estimate how many copies of the reporter are sufficient for quantifying the global transcription rate.To do so, we simulated the transcriptional dynamics of genes using two-state models in the presence of both gene-specific and gene-nonspecific (global) fluctuations.Gene-nonspecific fluctuations were implemented by assuming a fluctuating cellular biochemical environment that globally influences the rate of nascent RNA production for all genes (Materials and Methods).To illustrate the effect of both sources of fluctuations, we simulated the transcriptional dynamics of a gene under various global fluctuations, and found that the modulation by global fluctuations can be challenging to identify by looking at a single trace (Supplementary Figure S2B).By summing multiple reporter genes, we could successfully recover the globally modulated transcriptional activity (Figure 2 A bottom).We further found that ∼10 copies of the reporter gene would be necessary for capturing global transcription rate dynamics when the global fluctuations are in a sinusoidal waveform, and the fluctuation amplitude is ∼70% of the reporter expression le v el (Supplementary Figure S2C; Materials and Methods).We demonstrated that this result is generally robust to variations in the underlying kinetic model (e.g. when including two elongation reaction steps to account for the complexity in the RN A pol ymerase elongation dynamics ( 55-57 ) (Supplementary Figure S2D), or when increasing or decreasing the promoter switch rates (Supplementary Fig- ure S2E)) and the waveform of global fluctuations (Supplementary Figure S2F).

Multi-integr ated tr anscriptional r eporters allow ed f or temporally analyzing global transcription rate dynamics
Based on the in-silico results, we next aimed to construct and test the multi-reporter system experimentally.Our system is based on a synthetic transcriptional reporter containing CFP and RNA binding stem-loops (24 × PP7) ( 54 ), dri v en by doxy cy cline-inducib le promoter TRE3G.When transcribed, the 24 copies of the PP7 stem-loop would recruit stab ly e xpressed PCP-mCherry protein, allowing the visualization of nascent transcription sites as fluorescent puncta (Figure 2 A, top).This reporter was then multiintegrated into U2OS cells using PiggyBac transposase.The r esulting monoclones wer e scr eened by imaging, and a monoclone with ∼12 visible nascent transcription sites was chosen, which was then sequenced to map reporter integration sites.We identified nine integrated reporter copies across se v en different chromosomes by genome sequencing (Supplementary Figure S3A), yet there are likely more reporter copies that were not captured by sequencing (Supplementary Figure S3B).
Using the reporter cell line, we performed time-lapse imaging of single cells and quantified the fluorescence intensity traces of individual nascent transcription sites (i.e.individual fluorescent puncta in the nucleus).Intriguingl y, w hen summing transcriptional activity traces from all reporter loci in a single cell, w e w ere able to observe large fluctuations in the summed activity dynamics (Figure 2 B, Supplementary Figure S3C, Supplementary Video S1).This observation is consistent with the in-silico result that summing over a sufficient number of reporters can r ecover shar ed global fluctuations.Howe v er, such global fluctuations could result from fluctuations either inside the cell (e.g.changes in the intracellular environment) or outside the cell (e.g.changes in the culture condition).To distinguish between the two sour ces, we compar ed the summed transcriptional activity traces from cells in the same culture environment and found that global fluctuations were unsynchronized across cells, indica ting tha t the fluctua tions are cell-intrinsic (Supplementary Figure S3D).Together, these results indicated that the reporter system could capture global transcription rate dynamics in single cells.It should be noted that the summed transcriptional activity dynamics only provide an estimate of the global transcription rate dynamics, i.e. such a reporter system focuses on capturing the temporal pattern and the relati v e magnitude, instead of the absolute magnitude, of the global transcription rate.
Ne v ertheless, the abov e data prompted se v eral questions regarding the characteristics of global transcription rate dynamics.First, because the reporter genes are dri v en by the same promoter (i.e.doxy cy cline-inducib le promoter, TRE3G), it is critical to discern whether the observed dynamics are indeed global, instead of regulon-specific.Second, because the transcriptional dynamics of individual reporter genes are subjected to both gene-specific and gene-nonspecific fluctuations, it is necessary to compare the relati v e strength of the two noise sources inside cells.Additionally, because the magnitude of fluctuations should depend on the le v el of transcription, it is necessary to quantify such dependency.Third, besides characterizing the magnitude of fluctuations in global transcription ra te d ynamics, it is also critical to profile their temporal characteristics.
First, to discern if the fluctuations in the summed activity dynamics are global or regulon-specific, we needed to test whether genes dri v en by dif ferent promoters (i.e.dif ferent regulons) would be influenced by the same fluctuations.To do so, we constructed a cell line with two different types of transcriptional reporters (by using PP7 and MS2 systems ( 52-54 )), each dri v en by different promoters (e.g.TRE3G and AT CB promoter, respecti v ely) (Figure 2 C , Ma terials and Methods).If the fluctuations are regulon-specific, the dynamics of the two types of reporters would not exhibit a significant correlation, and vice versa if the fluctuations are global.As a negati v e contr ol, reporters fr om different cells should not correlate.By performing two-color time-lapse imaging, we successfull y ca ptured the transcriptional dynamics of two different types of reporters using two separate RNA-binding protein systems (Supplementary Figure S3E, Figure 2 D).We found that reporters from two regulons in the same cell exhibited a significantly positi v e correlation.In contrast, reporters from different cells were not correlated (Figure 2 E).This result is consistent with the presence of global fluctuations influencing regulons nonspecifically.Characterization of another cell line containing reporters dri v en by TRE3G and UBC pr omoters pr ovided additional support for this conclusion (Figure 2 E).
To further address whether it is sufficient to capture global fluctuations using reporter genes of a single regulon or if reporter genes from different regulons are necessary to remove potential 'regulon-specific' effect, we computationally compared differ ent r eporter designs, with reporter genes from one to three regulons (Supplementary Figure S3F).We found that the capability to capture the global fluctuations depended on the number of reporter genes, rather than the number of regulon types used (Supplementary Figure S3G).Importantly, for the design with two types of regulons (which is analogous to the tworegulon e xperiments abov e), the correlation between the summed transcriptional activities of the two regulons was markedly lower than the correlation between each regulon's summed transcriptional activities and the global fluctuation waveform (Supplementary Figure S3G), consistent with the relati v ely low correlation between regulons observed in the experiments (Figure 2 E).
We next addressed how such global fluctuations are compared to gene-specific (i.e.intrinsic) fluctuations in terms of the magnitude of fluctuations.To do so, we decomposed fluctuations in transcriptional activity into gene-specific and gene-nonspecific components by adapting the concept from the two-color reporter assay ( 2 , 42 ).More specifically, based on activity traces from any two reporter loci in a cell, we used mean-normalized variance and covariance between two signals to estimate gene-specific and genenonspecific sources of fluctuations, respecti v ely (Supplementary Figure S4A, Materials and Methods).By performing such calculations across all pairs of reporter loci in a cell and using the averaged results as cell-specific fluctuation strength, we found that the strengths of both fluctuations decrease as transcriptional activity increases (i.e.doxycycline concentration increases) and that gene-specific fluctuations are more substantial than gene-nonspecific fluctuations across different doxy cy cline conditions (Supplementary Figure S4B).Thus, genes transcribing at lower activity le v els ar e mor e susceptible to global fluctuations compared to those transcribing at higher activity levels.Moreover, because gene-specific fluctuations are much stronger than gene-nonspecific fluctuations, these results implicate that gene-specific control of transcriptional bursting (i.e. the source of gene-specific fluctuations) mainly contributes to the strength of stochasticity in individual genes' transcriptional activities, whereas gene-nonspecific fluctuations mainly contribute to the temporal organization of transcriptional activities (Supplementary Figure S4C), highlighting the importance of temporal analysis for studying gene-nonspecific fluctuations.
We further investigated the temporal characteristics of global tr anscription r a te d ynamics by performing autocorrelation anal ysis, w hich detects potential periodicity.We computed auto-correlation functions of individual cells' summed transcriptional activity traces under various doxycycline concentrations, and estimated the periodicity by measuring the distance between the first two peaks in the auto-correlation function.The results showed that global fluctuations consistently exhibit a ∼12 h periodicity across different transcriptional acti vity le v els (Supplementary Figure S4D, E), indicating that global fluctuations are potentially dri v en by or coupled to periodic cellular signals.
To gether, using the m ulti-integrated reporter system, we re v ealed that the gene-nonspecific fluctuations affecting reporter transcription are global instead of regulon-specific, and that such fluctuations are much weaker than genespecific fluctuations inside a cell.Importantly, the reporter systems allowed us to capture global transcription rate dynamics by averaging out gene-specific fluctuations, revealing the potentially periodic nature of such dynamics.

Global tr anscription r ate contains both static and dynamic components of non-genetic heterogeneity
While global transcription rate fluctuates temporally within a single cell, the mode of the resulting non-genetic heterogeneity within a clonal cell population remains to be further examined.This would allow us to address whether cells can access distinct or overlapping global transcription rates, and whether these rates are heritable between generations, which are necessary for further dissecting the mechanism and implications.
We first aimed to discern between static and dynamic modes of heterogeneity using dynamic data.More specificall y, we asked w hether the global transcription rates in individual cells are persistently and statically different (i.e.static heterogeneity) or whether different cells are dynamically accessing a shared distribution of ra tes (i.e.d ynamic heterogeneity).By examining the global transcription rate d ynamics of dif fer ent cells (Figur e 3 A), we found that while global transcription rate fluctuates within each cell, the accessible rates of different cells come from overlapping distributions (Figure 3 B).Moreover, the mean of each cell's rate (i.e.time-aver aged global tr anscription r ate) is much more widely distributed when compared to the scenario where all gene loci were independently and randomly regulated (Figure 3 C, Supplementary Figure S5A, Materials and Methods).These results indicate the co-existence of both dynamic and static hetero geneity, w hereby global transcription rates in individual cells fluctuate around distinct mean rates but with overlapping instantaneous rates.
The above findings prompted us to further characterize the dynamic aspect of the heterogeneity by depicting how the global transcription rate temporally fluctuates in individual cells.Because the global fluctuations exhibit a ∼12h periodicity (Supplementary Figure S4D), which is close to half of the cell cycle length, we speculated that global transcription rate could be dynamically coupled to the cell cycle (Figure 3 D).To explore this, we constructed a cell cycle sensor based on PCNA fusion protein (Citrine-PCNA) ( 58 ), allowing us to identify G1, S and G2 / M phases based on the Citrine-PCNA puncta formed during DNA replication (Figure 3 E).By imaging reporter transcription and citrine-PCN A signals sim ultaneousl y in the same single cells, we found that global transcription rate is indeed cell-cyclecoupled, eliciting two pulses in a cell cycle (Figure 3 F, Supplementary Figure S5B).More specifically, the rate le v el peaked at the G1 phase and late S or early G2 / M phase, respecti v ely, and reached the trough at the mid-S phase.Note that when reporters transcribed strongly (i.e.under the highest doxy cy cline concentration), the two pulses are not as obvious (Figure 3 F), consistent with the earlier finding that genes with higher transcriptional activities are less susceptible to global fluctuations (Supplementary Figure S4B).
To further dissect the coupling between the cell cycle and global tr anscription r ate, we asked whether perturbing the cell cycle would affect global transcription rate dynamics (Supplementary Figure S5C).We thus used thymidine to block DNA replication and found that the rate dynamics no longer exhibited a pulsatile behavior, but instead sustained at a relati v ely high le v el when cells were blocked at the S phase (Supplementary Figure S5D, E).This result suggested that the troughing of global tr anscription r ate is coupled to S phase progression, and that the global transcription rate is acti v ely modulated during the cell cycle.
Gi v en the co-existence of both dynamic and static heterogeneity, we next focused on the static aspect of the heterogeneity.Note tha t sta tic heterogeneity is used her e to r efer to variability between cells, while dynamic heterogeneity is used to refer to temporal variability within specific cells.While we have demonstrated that individual cells can ha ve distinct time-a veraged global transcription rates (Figure 3 C, Supplementary Figure S5A), indicati v e of static heterogeneity, it remains to be addressed whether such static heterogeneity is heritable across generations.We thus computed the correlation between cell-cy cle-av eraged global tr anscription r ates of mother and daughter cells, and observed a significant positive correlation (Figure 3 G).Note that 3-B5-6 / citrine-PCNA cell line (Materials and Methods) was used in this experiment (see Note S2 for more details).In other words, the global transcription rates in both mother and daughter cells would temporally fluctuate around mean le v els that are significantly correlated.Thus, the global transcription rate appears to be heritable between successi v e generations, which would contribute to the persistence of clonal heterogeneity.It should be noted that such heritability is likely transient, i.e.only across neighboring generations, as the R 2 value is not very high.
To further quantify heterogeneity, we next compared the contributions from dynamic and static aspects of heterogeneity to non-genetic heterogeneity in global transcription rate.Non-genetic heterogeneity is typically measured in a static population of cells (e.g. by flow cytometry, image snapshots, or scRNA-seq), which can be regarded as the total heterogeneity in a population and can be quantified by the coefficient of variation (CV).Our temporal measurements allowed deconvolving the total heterogeneity into ones contributed by different aspects, including dynamic and static heterogeneity (Materials and Methods).To quantify the degree of static heterogeneity, we computed the coefficient of variation (CV) of the population-averaged and cell-cycle-normalized global tr anscription r ate tr ajectory (Figure 3 H, Materials and Methods).To quantify the contribution of dynamic heterogeneity, we computed the CV of the time-averaged rates from all individual cells.We found that compared to dynamic heterogeneity, static heterogeneity yields a larger contribution to the total heterogeneity in global transcription rate (Figure 3 H), and that the relati v ely large residual CVs indicated the presence of fluctuations intrinsic to global transcription rate (that were not averaged away by the m ulti-reporter assay).Reassuringl y, simulations of a population of cells with intrinsically different mean global transcription rates supported key findings of the experimental results (Supplementary Figure S5F-I): (i) both static and dynamic heterogeneities can be consistentl y extracted w hen using differ ent r eporter designs (i.e.either the strength or the copy number of the reporter); (ii) the residual CV primarily represents the intrinsic noise of the estimated global transcription rate as it could be largely reduced by increasing the copy number.
Because cell size has been implicated in global transcription control ( 34 , 35 ), we next asked whether the observed heterogeneity in global transcription rate relates to cell size.
To do so, we used the size of the nucleus as a correlate of cell size ( 59 ) (Supplementary Figure S6A) and characterized how time-averaged nuclear size correlates with timeaver aged global tr anscription r a te, and how nuclear size d ynamically changes over the cell cycle in individual cells.Interestingly, we found that time-averaged nuclear size displays a statistically significant (but a relati v ely low magnitude) correlation with time-averaged global transcription rate across single cells onl y w hen the doxy cy cline concentration is low (Supplementary Figure S6B), indicating that nuclear size differences could contribute to the static heterogeneity of the global tr anscription r ate for lowly transcribing genes.Meanwhile, we found that the dynamic patterns of nuclear size along the cell cycle are distinct from those of global transcription rate (Supplementary Figure S6C).Ther efor e, further temporal studies are necessary to substantiate the role of cell size in global transcription.

Single-cell transcriptome analysis supported the presence of a hybrid mode of heterogeneity in the global transcription rate
Thus far, we have established that global transcription rate is heterogeneous in a clonal population of cells, and such heterogeneity is contributed both by static and dynamic components.Because these results were obtained using temporal data from the reporter system, we next sought to employ single-cell transcriptome data of the same cell type to validate key conclusions.
We ther efor e r esorted to the earlier da taset tha t contains transcriptome measurements of cell-cycle-sorted U2OS cells ( 45 ).This dataset would allow us to test whether cells from different cell cycle phases would exhibit variable global tr anscription r a tes (i.e.d ynamic heterogeneity), and if there ar e inter cellular variabilities within cells from the same cell cycle phases (i.e.static heterogeneity) (Figures 4 A, B).
By using summed intronic read counts (which are spikedin normalized) as a proxy of global tr anscription r ate and by binning cells based on cell cycle phases, we calculated the averaged global transcription rates along the cell cycle (Materials and Methods).The resulting 'pseudo dynamics' of the global tr anscription r ate exhibited two activity pulses within one cell cycle (Figure 4 C), w hich reassuringl y agrees with the temporal dynamics from our reporter system (Figur e 3 F).Mor eover, when comparing genes transcribing at dif ferent ra tes across the transcriptome, we found that lowexpr ession genes wer e mor e susceptible to the modulation (Figure 4 C).These results supported the presence of dynamic heterogeneity and further substantiated our conclusion that gene transcriptional activity is indeed globally and temporally modulated in a gene-nonspecific manner.
To quantify static heterogeneity, we needed to compute intercellular variability among cells within the same cell cycle phase and compare such variability to variability arising from Poisson fluctuations.We thus placed cells into 50 bins along the cell cycle, and within each bin, the CV of the total intronic read counts was calculated (Figure 4 D).To calculate CV from Poisson fluctuations, we modeled the detection of each gene's introns as a Poisson process, and the CV was then estimated as the inverse of the squared root of total intron counts.By doing so, we found that intercellular variability is considerably above Poissonian variability To detect dynamic heter ogeneity, cells fr om the different cell cycle phases could be separately averaged to obtain the dynamics along the cell cycle.To detect static heterogeneity, variability in cells from the same cell cycle phase could be compared to the variability of a control distribution ( A ).To achie v e this, cell-cy cle-sorted single-cell RNA-seq data from Mahdessian, 2021 was used ( B ) (with spiked-in RNAs).Intronic read counts summing over genes of differential e xpression le v els in indi vidual cells were calculated, and results from cells along different cell cycle bins were averaged ( n = 33-79 cells in each bin) ( C ).The variability of total intron counts (of low expressing genes) of cells in each cell cycle bin was compared with the variability of a Poisson distribution of cells with the same mean count ( D ).Note that in (D), the CVs from Poisson distributions are between 0.003 and 0.005 (which appeared to be zero in the plot).
(Figure 4 D), supporting the existence of intercellular variability in the global transcription rate.
To gether, anal yses of single-cell transcriptome data provided additional evidence supporting the co-existence of dynamic and static heterogeneity in global transcription rate, highlighting the perplexing nature of such non-genetic heterogeneity.

Perturbation assays suggested the role of ATP metabolism in conferring the non-genetic heterogeneity of global transcription rate
To further decode such heterogeneity, we next focused on its mechanistic underpinnings.We began by le v eraging single-cell transcriptome data to search for potential regulators whose expression levels can explain the differences in global transcription rates among single cells.In doing so, we assumed that both inter-and intra-cellular variabilities in global tr anscription r ate are potentially dri v en by changes in specific regulators that lead to changes in the global cellular state.To detect such regulators, we used a gene regulatory network inference package based on ran-dom for est r egr ession (i.e.GEINE3 ( 60 )) to rank genes based on their explanatory powers for global transcription rate (Figure 5 A).Intriguingly, the top 100 genes are highly enriched for biological processes related to ATP metabolism (Figure 5 B), indicating that ATP metabolism may be playing a key role in modulating global transcription rate.
Following this lead, we asked whether ATP le v el correlates with global transcription rate in single cells.We used an ATP biosensor ( 61 ) to measure nuclear ATP le v els and subjected the same cells to a pulse of nucleoside analog 5-EU (5-ethynyl uridine) to measure instantaneous global transcription activities (Materials and Methods).The rationale to use 5-EU labeling for measuring global transcription activity is that the transiently incorporated 5-EU can r epr esent the instantaneous transcription activity of all genes in the cell and that the cell line containing the ATP biosensor is technically challenging to install the PP7-based reporter system (i.e. it would r equir e generating a monoclonal cell line).Importantly, we demonstrated that the amount of incorpora ted 5-EU correla ted with the summed PP7 signals in the reporter cell line (Supplementary Figure S7A, B pro viding a k ey support for using either 5-EU labeling or PP7-based reporter system for quantifying global transcription activity.Cells with ATP biosensor were stained to label the incorporated EU, after which we sim ultaneousl y quantified ATP le v els and global transcription activity (estimated by the amount of incorporated EU) in the same cells with two-color fluorescence imaging (Supplementary Figure S7C).Intriguingly, we found that the two signal le v els are significantly correlated (Supplementary Figure S7D), consistent with a hypothesis that metabolites produced during ATP metabolism (such as ATP) could globally modulate gene transcriptional activity nonspecifically, leading to heterogeneity in the global transcription rate.
To test this hypothesis, se v eral key issues needed to be addressed.First, it is important to validate whether the global transcriptional modulation is mediated via a diffusible mechanism (e.g. through metabolites) or a nondiffusible mechanism (e.g. through chromatin dynamics).Second, if ATP metabolism plays a key role, it is necessary to test whether perturbations to ATP metabolism would impact the dynamics and the le v el of global transcription rate.Third, because the global transcription rate is heritable between successi v e generations, it is critical to deciphering whether ATP metabolic state is also heritable.
We first tested whether a diffusible mechanism could be responsible for global transcriptional modulation.We reasoned that if metabolites from ATP metabolic reactions (that occur in cytoplasm or mitochondria) diffuse into the nucleus to modulate transcription globally, two nuclei embedded inside the same cytoplasm should exhibit correlated global tr anscription r ate le v els and coor dinated global transcription rate dynamics (Figure 5 C).By using PEG to fuse the membranes of two cells and imaging the resulting fusion cells to quantify global transcription ra te d ynamics (Figures 5 D, E), we reassuringly found that the time-averaged global tr anscription r ates ar e significantly corr elated between two nuclei in individual fusion cells (Figure 5 F), and that their dynamics are also highly coordinated (Figure 5 G), supporting that a diffusible mechanism potentially mediates the global modulation of gene transcription in the nucleus.
We next perturbed ATP metabolism and quantified the impacts on the global transcription rate.Because ATP metabolism involves both gl ycol ysis and oxidative phosphorylation, we reasoned that inhibiting either stage of ATP metabolism should lower the stead y-sta te le v el of global tr anscription r ate (Figure 5 H).When a ppl ying 2-DG and oligomycin to separately inhibit glycolysis and oxidati v e phosphorylation and imaging cells under perturbation for at least one cell cycle, we found that time-averaged global tr anscription r ates wer e significantly r educed in either perturbation (Figure 5 I).Notably, global tr anscription r ate was more susceptible to gl ycol ysis inhibition, indicating a more important role of gl ycol ysis than oxidati v e phosphorylation in the modulation.
Moreover, because the global transcription rate is dynamically modulated during the cell cycle, we dynamically inhibited different stages of ATP metabolism along the cell cycle and quantified the immediate responses, allowing us to establish the dynamic role of ATP metabolism on global transcription during the cell cycle (Figure 5 H).By doing so, we unexpectedly found distinct dynamic roles of glycol-ysis and oxidati v e phosphorylation on global transcription ra te modula tion (Figure 5 J).More specifically, inhibition of gl ycol ysis led to a significant reduction in global transcription rate, which is independent of the cell cycle (Figure 5 J left, Supplementary Figure S8A, B).In contrast, inhibition of oxidati v e phosphorylation had a cell-cycle-dependent effect, as it reduced global transcription rate only after the S phase (Figure 5 J right, Supplementary Figure S8C, D).These intriguing perturbation outcomes are consistent with a picture that the global tr anscription r ate switches to an increased dependency on oxidati v e phosphorylation for ATP metabolism during the mid-S phase while maintaining a steady dependence on gl ycol ysis throughout all cell cycle phases.
Based on these observations, we hypothesized that blocking cells at S phase would sensitize the response to the inhibition of oxidati v e phosphorylation, but not to the inhibition of gl ycol ysis.We thus used thymidine to block cells prior to the addition of inhibitors, and found that when blocked at the S phase, global transcription rates were indeed more sensiti v e to the addition of oligomy cin (Supplementary Figur e S8E-G) compar ed to the scenario without blocking (Figure 5 J right, Supplementary Figure S8C, D).In contrast, the responses to the inhibitor of gl ycol ysis did not display such a contrast (compare Supplementary Figure S8H-J with Figure 5 J, left and Supplementary Figure S8A, B).Such a differential sensitivity to gl ycol ysis inhibition versus oxidative phosphorylation inhibition can be further visualized when plotting all the results together (Supplementary Figure S8K).
Thus, these d ynamic perturba tion results implica ted tha t a cell-cycle-coupled switch in ATP metabolism (i.e.switch to an increased dependency on oxidati v e phosphorylation) modulates global transcription rate dynamics, whereby the decreasing energy supply (or other metabolites from ATP metabolism) prior to the switch dri v es the decrease in global tr anscription r ate, which is then boosted upon switching (giving rise to the second pulse of global transcription rate).
Lastly, a long-time imaging experiment using the ATP biosensor was performed, and the result showed that the time-av eraged ATP le v el was heritab le across successi v e generations (Supplementary Figure S8L), supporting the potential role of ATP metabolism in conferring heterogeneity in global transcription rate.
In sum, these results suggest an intriguing role of ATP metabolism in both the temporal modulation of global transcription rate as well as the modulation of time-averaged global transcription rate inside single cells, highlighting the potential function of metabolic dynamics in conferring non-genetic heterogeneity in mammalian cells.Howe v er, it should be noted that drug perturbations to metabolism could cause alerted metabolism of other metabolites, and we cannot exclude the roles of such alterations in global transcription rate.

Heterogeneous global transcription rate is associated with the intrinsic drug tolerance of cancer cells
Non-genetic heterogeneity of cancer cells has attracted growing interest, as it has been implicated in chemotherapeutic resistance , metastasis , and tumor microevolution ( 12 , 16 ).While global transcriptional amplification has been observed in cancer cells, most studies focused on the role of induced global transcriptional amplification ( 25 , 26 , 29 ), rather than the pre-existing heterogeneity in the global transcription rate.Nonetheless, as we have established that the global transcription rate is heterogenous at steady state (without perturbation) and is transiently heritable across cell generations, it remains elusi v e whether such pre-existing intrinsic non-genetic heterogeneity may have functional implications.
Because of the earlier study indicating the potential role of global transcriptional amplification in drug resistance ( 29 ) as well as the r ecent inter ests in the role of pre-existing non-genetic heterogeneity in chemotherapeutic resistance ( 30 , 62 , 63 ), we set out to characterize the functional implication of heterogeneous global transcription rates in intrinsic drug tolerance.We reasoned that because individual cells have variable time-aver aged tr anscription r ates prior to drug treatment, such variability could relate to differential drug tolerances for different cells (Figure 6 A).To test this, we carried out a single-cell drug tolerance assay to quantify the relationship between drug tolerance and global transcription rate (Supplementary Figure S9A).
More specifically, we imaged cells with the reporter system for 24 h prior to the addition of a commonly used chemotherapeutic drug, etoposide (targeting topoisomerase II), and continued the imaging for another two days (Supplementary Figure S9B, Materials and Methods).By quantifying the intrinsic global transcription rates in individual cells and measuring the survival time of the corresponding cells in the presence of etoposide, we found that the intrinsic global transcription rate exhibited a statistically significant positi v e correlation with the survival time (Supplementary Figure S9C), yet the magnitude of the correlation is relati v ely moderate, possib ly because of both the biological noise in the system as well as the technical challenge in accurately determining the time of cell death.Intriguingly, similar r esults wer e observed for two additional chemotherapeutic drugs, bortezomib (inhibiting the 26S proteasome) (Supplementary Figure S9D) and cisplatin (DNA damaging) (Supplementary Figure S9E).
Due to the limited throughput of imaging-based survival assay, we next aimed to test more drugs by developing an assay with higher throughput.We enriched cells with low or high global transcription rates by FACS sorting after ∼12 h of doxy cy cline induction (Supplementary Figure S9F, Materials and Methods) and then compared their susceptibility to drugs.Reassuringly, this new assay recapitulated results from the imaging assay, i.e. cells with higher global transcription rates displayed higher tolerance to the three drugs tested above (the leftmost three drugs in Figure 6 B, Supplementary Figure S9G).Note that 3-B5-6 / Citrine-PCNA cell line (Materials and Methods) was used in this experiment (see Note S2 for more details).We then used this assay to study 26 more drugs, and among the total of 29 drugs, only 17 drugs display appropriate cell-killing effects (Materials and Methods, Supplementary Table S1).Among these 17 drugs, six were significantly better tolerated by cells with higher global tr anscription r ates (Figure 6 B), and the rest (except one) were not dif ferentially tolera ted (Supplementary Figure S9H).These results are generally consistent with the picture that cancer cells with high global transcription rates exhibited increased tolerance to certain anti-tumor drugs, yet the difference among drugs remains to be explored.While this finding reconciles with the recent report that drug-induced (i.e.non-intrinsic) global transcriptional amplification is relevant to drug resistance ( 29 ), our result highlights the potential implication of pre-existing (i.e.intrinsic) heterogeneous global transcription rates in drug tolerance.Note that these results were based on the 3-B5-6 / Citrine-PCNA cell line, which was constructed by lentiviral transduction of a cell cycle reporter (i.e.Citrine-PCNA) into a monoclonal transcriptional reporter cell line (i.e.3-B5-6).While we cannot rule out the potential caveats arising from the vir al tr ansduction step, we belie v e it is reasonab le to assume the monoclonality of the transcriptional reporters, as Citrine-PCNA introduced to the 3-B5-6 monoclone is irrelevant to the transcriptional reporter system (see Note S2 for more details).
To further explore this implication, we resorted to cancer cell line databases, including CCLE and GDSC (64)(65)(66)(67), where the bulk transcriptome profiles and drug sensiti vity profiles hav e been characterized for an e xtensi v e collection of cell lines and drugs.The rationale is that different cell lines may possess differential global transcription rates and thus variable intrinsic drug sensitivities (Supplementary Figure S10A, B).We sought to deri v e a marker gene set whose expression level could be used to estimate global tr anscription r ates in specific cell lines.Using the top 100 genes identified in our previous analysis (that can best explain global transcription rate in the single-cell sequencing dataset) (Figure 5 B) and the expression profiles of the cancer cell lines from the database, we identified the largest coexpressing gene module consisting of 28 genes.We used this gene set to compute a modulation index for each cell line as a proxy for the extent of global transcriptional modulation (Supplementary Figure S10C, Materials and Methods).
We first investigated drugs in the da tabase tha t overlap with the 17 drugs that were experimentally measured.Reassuringly, for the three overlapping drugs, the drug sensitivity data from the database agree with our experiments (Figur e 6 C).Mor e specifically, for two of the overlapping drugs (bortezomib and cisplatin), cell lines with high modulation indexes displayed higher tolerances compared to cell lines with low modulation indexes, consistent with our experimental observations (Figure 6 C, B, Supplementary Figure S9C-E); and for the third overlapping drug (OF-1), no differential drug tolerance was observed for both the database and our experiment (Figure 6 C, Supplementary Figure S9H).
We next asked whether drug sensitivity would display an overall correlation with modulation index in the database.Using the drug sensitivity data for camptothecin (a topoisomerase inhibitor) as an example, we found that the modulation indexes in different cell lines exhibit a significant positi v e correlation with either of the two drug sensitivity measures, IC50 and AUC (Supplementary Figure S10D, E).Importantly, similar r esults wer e observed for most of the drugs targeting different signaling mechanisms (Supplementary Figure S10F, G), indicating that the modulation index of a cancer cell line relates to the degree of  S1 for information on concentrations).n = 3-10.P -value from Welch's t -test.( C ) Bo xplots sho wing AUC values (characterizing drug tolerance le v el) from GDSC2 database for cell lines with low or high modulation indexes.n = 39 and 30 cell lines with low or high modulation indexes, respecti v ely.P -value from Welch's t -test.intrinsic multi-drug tolerance of the cell line.Ther efor e, this result further indicates the potential biological implication of global transcription rate modulation.It should be noted that the above analyses using the public cancer cell line databases were built on the assumption that the modulation index correlates with the extent of global transcriptional modulation, and thus the implications based on such analyses would need further investigations.

DISCUSSION
It has been well established that the expression levels of individual genes are often highly heterogeneous in a clonal population, and such heterogeneity can play pivotal roles in a range of biological processes.In contrast, due to technical challenges, cell-to-cell variability in the expression of the entire transcriptome has remained largely unexplor ed.We addr essed this challenge by integrating timelapse imaging with transcriptome analysis, and uncovered that gene transcription is globally and temporally modulated in a gene-nonspecific manner, resulting in heterogeneous and heritable global transcription rates in individual cells (Figure 7 ).Intriguingly, our single-cell drug tolerance assay and analyses of cancer cell line datasets indicated the potential implication of such heterogeneity in conferring intrinsic drug tolerance, underscoring the importance of studying non-genetic heterogeneity at the global scale.Deciphering the sources of fluctuations is critical for understanding heterogeneity in a clonal population.Because global (gene-nonspecific) fluctuations can convolute with gene-specific fluctuations to confer heterogeneity in a clonal population, it has thus been highly challenging to quantify the global source of fluctuations ( 18 ).Our reporter system was rationally designed to le v erage the correlated nature of global fluctuations among genes, allowing us to temporally analyze fluctuations in global transcription ra te tha t af fect all genes nonspecifically (Figure 2 ).Compared to temporal analysis using our reporter system, although single-cell transcriptome analysis can only provide a snapshot of cells and thus cannot distinguish between modes of heterogeneity, it can still be used to infer the presence of transcriptome size heterogeneity ( 43 , 44 ).Thus, as a complement to our temporal analysis, it would be helpful to systematically quantify absolute single-cell transcriptomes across di v erse conditions and cell types to determine the pervasi v eness of heterogeneity in the global transcription rate.Howe v er, it is important to note that most existing single-cell transcriptome datasets cannot be used for analyzing such heterogeneity, as protocols without using external spiked-ins do not measure absolute transcriptomes.
While genes are often specifically regulated by upstream regulators to cope with various environments, increasing evidence has suggested that the activities of many genes can also be globally modulated during important cellular processes, such as cancer progression and stem cell differentiation ( 68 , 69 ).Our finding provides quantitati v e e vidence that global modulation is not only highly dynamic, but also highly cell-specific, displaying features of both static and dynamic heterogeneity (Figure 3 and Supplementary Figure S5).Notably, mitochondrial ATP production has been shown to be a global influencer of gene transcription ( 37 ), and it was thought that the cell-to-cell variability in global gene transcription is contributed by the variability in the mitochondrial mass ( 37 ).Howe v er, while mitochondrial mass variability may provide an explanation for the static component of heterogeneity we observed in the global transcription rate, se v eral lines of evidence from our study point to the potential role of dynamic ATP metabolism in shaping such heterogeneity (Figure 5 and Supplementary Figure S8), underscoring the emergent roles of cellular metabolism in the temporal organization of cellular transcriptome.Yet, because most of the key evidence was deri v ed from perturbation experiments using 2-DG or oligomycin, we could not exclude the potential roles of other metabolites perturbed in the experiments.Ther efor e, how ATP metabolism dynamically modulates global tr anscription r ate and what specific and fundamental molecular mechanisms are involved in the modulation r equir e further investigations, so as the role of additional factors or processes other than ATP metabolism.More generally, the interplay between metabolism, cell size, and global transcription ra te modula tion a t the single-cell le v el necessitates further studies.
At the functional le v el, dynamic and cell-specific modulation of the transcriptome may bias cells towards specific states or fates while maintaining dynamic plasticity.As we have demonstrated, cells with high global transcription rates were able to survi v e longer when challenged with three different chemotherapeutic drugs (Figure 6 ).While this observed correlation implicates the potential role of heterogeneous global transcription rate and is consistent with the result from the cell line data analysis, it remains elusi v e what molecular entities are responsible for the longer drug survival and what other functional roles the heterogeneity may play.Nonetheless, our result highlights the necessity to carefully examine the consequences of variability within a clonal population, and adds to the growing list of evidence supporting the important, but often hidden, role of non-genetic heterogeneity in cancer cells ( 12 , 24 , 30 , 62 ).

DA T A A V AILABILITY
All relevant data have been included in this manuscript.Additional data are available from the corresponding author upon request.

SUPPLEMENT ARY DA T A
Supplementary Data are available at NAR Online.

Figur e 1 .
Figur e 1. Non-genetic hetero geneity in the global transcription ra te.(A, B) Quantifica tion of cell-to-cell variability in transcriptome size using spiked-in normalized single-cell RNA-seq data.( A ) Distribution of spiked-in normalized total exon counts in single U2OS cells ( n = 986 cells).Single-cell RNA-seq data from Mahdessian ( 45 ) was used (see Materials and Methods for details).( B ) Scatter plots showing total exon count per cell versus total intron count per cell (left) or the ratio of total intron count over total exon count (right).Total intron count was used as a proxy of the rate of global transcription, while intron / exon ratio is a proxy for global mRNA degradation rate.See Supplementary Figure S1A for explanations.( C ) Analyses of snapshot single-cell RNA-seq data pro vok ed questions regarding the mode, mechanism, and implication of cell-to-cell variability in the global transcription rate in a clonal population (i.e.non-genetic heterogeneity).

Figure 2 .
Figure 2. De v eloping a reporter system for the temporal analysis of global tr anscription r a te.( A ) Design schema tic (top) and in silico valida tion of the multi-integrated reporter system (bottom).See Materials and Methods and Supplementary Figure S2 for details on sim ulations.( B ) Time-la pse imaging of the reporter cell line allows temporal analysis of global transcription rate.An example cell labeled with identified nascent transcriptional sites (left), tempor al tr anscriptional activity dynamics of these sites (middle), and two dif ferent quantifica tions of global transcription ra te d ynamics (right) were shown.Example snapshots of the example cells were shown in Supplementary Figure S3C.See also Supplementary Video S1. (C-E) Assay for testing whether the observed fluctuations are regulon-specific or shared between regulons.In this assay, two ortho gonal RN A-binding protein-based transcriptional reporter systems were used to report the transcriptional dynamics of two different regulons in the same single cells ( C ). Transcriptional activity dynamics of two separate regulons (i.e.PP7-based reporters dri v en by TetO promoter and MS2-based reporters dri v en by the promoter of ACTB gene) were quantified in individual cells by summing over the transcriptional activities of detected loci of each regulon ( D ).Correlation between transcriptional dynamics of differ ent r egulons in the same single cells was significantly higher than the correlation computed using scrambled cell identities ( E ).Two separate cell lines (i.e.pTetO-PP7 / pACTB-MS2 and pTetO-PP7 / pUBC-MS2) were used for the quantification ( n = 51 and 33 cells, respecti v ely).P -value was from Welch's t -test.

Figure 3 .
Figure 3. Multi-integrated real-time transcriptional reporters allowed characterizing non-genetic heterogeneity in the global transcription rate.(A-C) Inferring the mode of heterogeneity using temporal data from the reporter system.Global tr anscriptional r a te d ynamics from three example single cells acquired under indicated condition ( A ) were used to construct distributions of instantaneous rate le v els ( B ).The distribution of time-averaged rate le v els in individual cells (n = 51 cells) was compared to a control distribution (i.e.binomial, see Materials and Methods) ( C ). (D-F) Characterizing the dynamic mode of heterogeneity in individual cells.To address whether global transcription rate is coupled to the cell cycle (left, D ), simultaneous quantification of the reporter system and the cell cycle is needed (right, D).Snapshots from simultaneous imaging of transcription and cell cycle marker PCNA-YFP at indicated time points post mitosis in an example cell ( E ).Population-averaged global transcription ra te d ynamics and the cell cycle signal (quantified by variance in nuclear PCNA signal) were plotted for fiv e indicated conditions ( F ). ( G ) Scatter plot showing time-averaged global transcription rate le v els in pairs of mother and daughter cells.( H ) Clonal variability in global transcription rate le v el from all cells in snapshots (characterized by the coefficient of variation, or total CV; 1st panel) can be decomposed into dynamic variability along the cell cycle time (characterized by dynamic CV; 2nd panel), static or inter-cellular variability (characterized by static CV, 3rd panel), and the residual CV (last panel).See Materials and Methods for details.

Figure 4 .
Figure 4. Single-cell transcriptome data supported a hybrid mode of heterogeneity in the global transcription rate.(A-D) Quantifying non-genetic heterogeneity in global transcription rate using spiked-in normalized single-cell transcriptome data (public data).To detect dynamic heter ogeneity, cells fr om the different cell cycle phases could be separately averaged to obtain the dynamics along the cell cycle.To detect static heterogeneity, variability in cells from the same cell cycle phase could be compared to the variability of a control distribution ( A ).To achie v e this, cell-cy cle-sorted single-cell RNA-seq data from Mahdessian, 2021 was used ( B ) (with spiked-in RNAs).Intronic read counts summing over genes of differential e xpression le v els in indi vidual cells were calculated, and results from cells along different cell cycle bins were averaged ( n = 33-79 cells in each bin) ( C ).The variability of total intron counts (of low expressing genes) of cells in each cell cycle bin was compared with the variability of a Poisson distribution of cells with the same mean count ( D ).Note that in (D), the CVs from Poisson distributions are between 0.003 and 0.005 (which appeared to be zero in the plot).

Figure 5 .
Figure 5. Evidence indicating the role of ATP metabolism in global tr anscription r a te modula tion.(A, B) Inferring potential regula tors of the global tr anscription r ate using single-cell transcriptome data.Single-cell expression matrix (normalized exon counts of genes by cells) was used as an input of GEINE3 for the inference of regulators for the global transcription rate ( A ). Results from gene ontology analysis of the top 100 potential regulator genes ( B ). (C-G) Cell fusion assay for testing if gene transcription is globally modulated by a diffusible factor.Design of the cell fusion assay ( C ). Snapshots showing two PEG-fused nuclei at indicated time points ( D ).Heatmaps showing z-scaled global transcription rate dynamics of two fused nuclei ( E ).Scatter plots showing two nuclei' time-aver aged r ate le v els in each single cell ( F ). Bo x plots sho wing the correlations between the dynamics of fused nuclei pairs and between the dynamics of non-fused nuclei pairs (i.e.control) ( n = 30 nuclei pairs and P -value from Welch's t -test) ( G ). (H-J) Responses of global tr anscription r a te to perturba tions in ATP metabolism.Assay includes two different experiments to measure the stead y-sta te responses and transient responses to inhibition of gl ycol ysis (using 10 mM 2-DG) or inhibition of oxidati v e phosphorylation (5 M oligomycin) ( H ). Boxplots showing single-cell time-aver aged global tr anscription r ate le v els under three conditions ( n = 36, 20, 23 cells from left to right, P -values from Welch's t -test) ( I ).Heat maps showing z-scaled global transcription ra te d ynamics in individual cells before and after indicated perturbations ( n = 27, 49 cells from left to right) ( J ).Cells were vertically aligned according to their cell cycle phases at the time of perturbation.See also Supplementary Figure S8A-D for additional quantifications.

Figure 6 .
Figure 6.The functional implication of intra-clonal global transcription rate variability in intrinsic drug tolerance.( A ) Hypothetical cartoon.Heterogeneous and heritable time-averaged global transcription rates in single cells could confer intrinsic differences in drug tolerance.( B ) Boxplots showing cell fraction survi v ed after indica ted drug trea tments for U2OS cells with low or high CFP le v els.Cells were sorted based on CFP intensity, subjected to drug treatments, and then quantified by imaging (see Materials and Methods for details and Supplementary TableS1for information on concentrations).n = 3-10.P -value from Welch's t -test.( C ) Bo xplots sho wing AUC values (characterizing drug tolerance le v el) from GDSC2 database for cell lines with low or high modulation indexes.n = 39 and 30 cell lines with low or high modulation indexes, respecti v ely.P -value from Welch's t -test.

Figure 7 .
Figure 7.A global source of non-genetic heterogeneity: mode, mechanism, and implication.Cartoons illustrating a working model for a global source of non-genetic heterogeneity in mammalian cells.Our results showed that ATP metabolism likely plays a key role in conferring heterogeneity in the global tr anscription r ate (left).In this model, cell cycle phase-dependent alter a tion in oxida ti v e phosphorylation acti vity may result in dynamic global transcription rates along the cell cycle, with the time-aver aged r ates being modulated by mean le v els of ATP metabolism.Intra-clonal variability in global transcription rate thus displays features of both dynamic and static modes of heterogeneity, with the static heterogeneity being heritable across successi v e generations (middle).Notab ly, such pre-e xisting heterogeneity between clonal cells could lead to differ ential r esponses to perturbations, and may thus have implications in intrinsic drug tolerance, among others (right).