Skip to main content
Download PDF

Fluorescence activated cell sorting followed by small RNA sequencing reveals stable microRNA expression during cell cycle progression

BMC Genomics volume 17, Article number: 412 (2016) Cite this article

Abstract

Background

Previously, drug-based synchronization procedures were used for characterizing the cell cycle dependent transcriptional program. However, these synchronization methods result in growth imbalance and alteration of the cell cycle machinery. DNA content-based fluorescence activated cell sorting (FACS) is able to sort the different cell cycle phases without perturbing the cell cycle. MiRNAs are key transcriptional regulators of the cell cycle, however, their expression dynamics during cell cycle has not been explored.

Methods

Following an optimized FACS, a complex initiative of high throughput platforms (microarray, Taqman Low Density Array, small RNA sequencing) were performed to study gene and miRNA expression profiles of cell cycle sorted human cells originating from different tissues. Validation of high throughput data was performed using quantitative real time PCR. Protein expression was detected by Western blot. Complex statistics and pathway analysis were also applied.

Results

Beyond confirming the previously described cell cycle transcriptional program, cell cycle dependently expressed genes showed a higher expression independently from the cell cycle phase and a lower amplitude of dynamic changes in cancer cells as compared to untransformed fibroblasts. Contrary to mRNA changes, miRNA expression was stable throughout the cell cycle.

Conclusions

Cell cycle sorting is a synchronization-free method for the proper analysis of cell cycle dynamics. Altered dynamic expression of universal cell cycle genes in cancer cells reflects the transformed cell cycle machinery. Stable miRNA expression during cell cycle progression may suggest that dynamical miRNA-dependent regulation may be of less importance in short term regulations during the cell cycle.

Background

The fine-tuned mechanisms of cell cycle have always been in the focus of cancer research resulting in a better understanding and optimization of the action of several chemotherapeutic agents [1–3]. Both posttranslational modifications of proteins (e.g. protein-protein interactions, phosphorylations) and altered transcriptional activity of specific genes contribute to the tightly controlled regulation of the cell cycle [3]. Analysis of mRNA transcripts expressed in a cell cycle dependent manner using high-throughput screening methods have identified numerous genes used for differentiating malignant tumors from benign lesions [4–6]. These correlations were presumably due to the possibility that cell cycle dynamics accelerate in malignant tumors, resulting in a larger proportion of cells residing in S and G2 phases [4].

Former approaches for detecting transcripts expressed in a cell cycle dependent manner used several synchronization techniques for arresting the cell cycle at a certain point. Among others, serum starvation, double thymidine block, and thymidine-nocodazole block halted the cell cycle in cultured cells at G0, early S and M phases, respectively [4–6]. After removal of the synchronizing agent, time-course gene expression data followed by adequate bioinformatics analysis were used for identifying the cell cycle regulated transcripts [4–6]. However, several conflicting arguments have been raised concerning the usage of these synchronization procedures [7]. Statistical re-examination of a former study [6] surprisingly revealed that randomization of time-course gene expression data showed the same strong periodicity in the expression patterns as those obtained from the original synchonization experiment. Moreover, synchronization procedures in general [8], and DNA replication inhibitors such as thymidine in particular, result in the perturbance of cell cycle machinery, producing growth imbalance and unscheduled expression of cyclins [8, 9]. Additionally, cells may lose their synchronization relatively soon after the release from the synchronizing agent [5, 10] and only a subset of cells reenter the cell cycle after arrest [5, 11]. Therefore, a set of criteria has been introduced for the analysis of cell cycle dependent transcripts [7]. Accordingly, the expression pattern of a certain gene can be introduced as cell cycle dependent if (i) no inhibition or starvation method was used for synchronization, (ii) results are reproducible over several experiments (iii) results of additional methods other than microarrays (e.g. Northern blot, quantitative real-time polymerase chain reaction – qRT-PCR) support the findings (iv) expression patterns are confirmed in non-synchronized experiments (e.g. cells separated by size or DNA content) and (v) statistically robust analysis supports the results [7].

The regulation of the cell cycle in general and the cell cycle dependent transcriptional program in particular is the consequence of the precise interactions between cyclin-cyclin dependent kinase complexes and an oscillating network of transcription factors [12, 13]. Additionally, epigenetic mechanisms as microRNA-mediated regulations contribute to proper cell cycle regulation. MicroRNAs (miRNAs) are short, ~22 nt long noncoding RNA molecules regulating gene expression on the post-transcriptional level targeting the 3’ untranslated regions of mRNAs [14]. Extensive complementarity results in mRNA degradation, whereas in the case of short complementarity, transcriptional silencing is achieved by transcriptional repression [15]. Among other physiological functions, the importance of miRNA-dependent gene regulation has been confirmed in several key members of the cell cycle machinery [16–18], contributing to miRNA-dependent cell cycle changes [17, 19]. Altered expression of cell cycle-controlling miRNAs has been reported in neoplasms of various tissues [18, 20, 21]. Additionally, dynamic miRNA expression changes have been observed during exit from quiescent state due to serum reintegration into the culture medium of serum starved cells [22, 23]. In particular, elevated expression of E2F1 and E2F3 in response to mitogenic stimuli have been shown to enhance the expression of its transcriptional targets: hsa-let-7 and hsa-miR-16 family members [22]. Moreover, E2F1 has been shown to enhance hsa-miR-15 expression, which inhibits cyclin E, one of the key transcriptional targets of E2F1 [17]. Accordingly, it has been proposed that such feed-forward loops encompassing the E2F transcription factors, miRNAs and cyclins contribute to the fine-tuning of cell cycle regulation [17]. However, the potential dynamic miRNA expression changes between the cell cycle phases of actively cycling cells without any synchronization or serum shock procedures have not been thoroughly investigated.

Here, for the first time, we show that gene expression signature obtained from unperturbed cells sorted by fluorescence activated cell sorting (FACS) based on their DNA content at different phases of the cell cycle correlate well with former gene expression studies using synchronization methods. In addition to lower expression of cell cycle genes in different cell cycle phases, dynamic mRNA expression changes were found to be of greater amplitude in primary, untransformed fibroblasts as compared to those detected in cancer cell lines, reflecting the more precise cell cycle regulation in untransformed cells with lower proliferation characteristics. Using numerous high-throughput miRNA-screening methods, miRNA expression, unlike mRNA expression, was found to be quite stable throughout the cell cycle progression in various human cells.

Results

An optimized cell cycle sorting method successfully differentiates cell cycle phases in various cells

Our optimal cell cycle sorting was able to differentiate cells residing in various cell cycle phases in all of the three cell types used (HDFa, NCI-H295R and HeLa cells) (Fig. 1, panel a-c). The purity of cell cycle sorted populations varied between cell types and cell cycle phases (Additional file 1: Table S1), but based on FACS reanalysis, these sorted cell populations were still more homogenous than cells obtained after synchronization procedures [5]. G1 phase was sorted most efficiently in all cell types with more than 95 % of purity in all cells. In NCI-H295R and HeLa cells, S phase cells showed more homogenous population as compared to cells in G2 phase.

Fig. 1
figure 1

Fluorescence activated cell sorting (FACS) analysis of stained cells and validation of cell cycle sorted populations in HDFa, NCI-H295R and HeLa cells. Panels a-c: FACS analysis shows cell cycle distribution of stained cells before (first image) and after sorting to G1, S and G2 phases (second, third and fourth images, respectively; Panel a: HDFa, Panel b: NCI-H295R, Panel c: HeLa. Intervals of fluorescence intensity (shown as vertical bands) were defined to gate G1, S and G2 phases, respectively. Panels d-e: Western blot analysis confirms the high efficacy of cell cycle sort (Panel d: NCI-H295R, Panel e: HeLa). Density values of phospho(Tyr15)-CDC-2 were first normalized to β-actin loading control and were further normalized to G1 phase (displayed above each band)

Full size image

Optimization of the cell cycle sort was needed for the sake of achiveving high purity in sorting without damaging or perturbing physiological cell functions. In particular, the determination of upper limit of sorting time, the use of a specialized sort medium and the immediate re-analysis of sorted cells contributed to our results.

Protein expression change confirms the successful sorting

Tyr15 phosphorylation of CDC-2 protein is a tightly controlled event in cell cycle progression [24], thus the respective amounts of phospho (Tyr15)-CDC-2 provide a general hallmark for each phase. Western blot analysis performed on protein extracts of sorted NCI-H295R and HeLa cells showed the well known phosphorylation patterns of CDC-2 (Fig. 1, panel d-e), confirming the purity of cell cycle sorting on protein level as well.

Gene expression profiling detects cell cycle regulated transcripts in FACS-sorted cells

The quantity and quality of isolated RNA from sorted cells were sufficient to perform high throughput gene expression screening (Additional file 2: Figure S1, Additional file 1: Table S1). Gene expression profiling, followed by rigorous statistical analysis detected 55 mRNA transcripts in NCI-H295R cells (Fig. 2, Panel a, Additional file 1: Table S3, panel B) and 252 mRNA transcripts in HeLa cells (Fig. 2, Panel b, Additional file 1: Table S3, panel C) to be expressed in a cell cycle dependent manner. Note that the majority of detected gene expression changes share a common manner: expression rises as cell cycle proceeds. Additionally, clustering showed that S and G2 phases’ expressional patterns are closer to each other than to G1 phase. Statistical analysis of HDFa microarray data failed to detect genes with significantly altered expression, however, the functional bioinformatics analysis of the gene expression changes of greater than a twofold change (FC > 2, Fig. 2, Panel f, Additional file 1: Table S3, panel A) supports the concept that these changes strongly influence cell cycle progression. Moreover, the successful qRT-PCR validation of the chosen FC > 2 genes (Fig. 2, Panel c) and the significant correlation of gene expression changes of cell cycle sorted with former synchronization-based experiments in primary fibroblasts (Fig. 3, Panel b) further confirm the relevance of our approach.

Fig. 2
figure 2

High throughput gene expression profiling, qRT-PCR validation and functional bioinformatics analysis of cell cycle dependent transcription in HDFa, NCI-H295R and HeLa cells. Panels a-b: Heat map of mRNA transcripts with significantly different expression between cell cycle phases (Panel a: NCI-H295R, Panel b: HeLa). Above the heat maps, hierarchical clustering upon differently expressed genes are also shown (additional information, including the list of gene symbols and descriptions can be found in Additional file 1: Table S3, panels B and C). Panels c-e: qRT-PCR validation of six genes chosen upon microarray analysis (Panel c: HDFa, Panel d: NCI-H295R, Panel e: HeLa). In each case ΔCt value was normalized to G1 phase (ΔΔCt) and was subjected to FC = 2-ΔΔCt transformation. Error bars show standard deviation. Asterisks mark statistical significance (p < 0.05). Panels f-h: Molecular and cellular functions concerned by gene expression alterations in cell cycle phases. Δ(G2-G1) gene expression changes of significantly differently expressed genes (Panel g: NCI-H295R, Panel h: HeLa) or genes with fold change > 2 expression (Panel f: HDFa) were subjected to IPA core analysis. In each cell, the five most significantly concerned networks are shown. Significance threshold (p = 0.05) corresponds to –log(p-value) = 1.301. For additional molecular and cellular functions see Additional file 2: Figure S2

Full size image
Fig. 3
figure 3

Comparison of cell cycle dependent gene expression observed by cell cycle sort and synchronization experiments and analysis of gene expression dynamics during the cell cycle in various cell types. Panel a: Venn diagram of cell cycle dependent genes detected in cell cycle sorted primary fibroblasts (HDFa SORT), in synchronized primary fibroblasts (PF synchr – data from [5]), in cell cycle sorted HeLa cells (HeLa SORT) and in synchronized HeLa cells (HeLa synchr – data from [4]). Intersections present the number of commonly found genes. For gene lists see Additional 1: Table S4. Panels b-c: Correlation analysis of gene expression differences using normalized expression values obtained from microarray experiments. Pearson’s correlation coefficient was calculated from Δ(G2-G1) expression changes detected by cell cycle sort and former synchronization experiments in primary fibroblasts (synchronization method: serum starvation – SS, Panel b) and HeLa (synchronization method: double thymidine block – DT, Panel c) cells. Correlation coefficients are displayed. Asterisks mark statistical significance (p < 0.05). For additional correlation calculations see Additional 2: Figures S3 and S4. Panels d-e: Analysis of normalized gene expression values of different cell types. Normalized gene expression of 127 cycling genes exported from microarray data (Panel d) and 10 cycling genes exported from qRT-PCR data (Panel e) were analyzed. Note the lower ΔCt(ACTB) values mean higher expression values. Error bars show standard deviation. Asterisks mark statistical significance (p < 0.05). Panels f-g: Analysis of mean gene expression fold change between cell cycle phases of different cell types. Absolute values of fold change of 127 cycling genes exported from microarray data (Panel f) and 10 cycling genes exported from qRT-PCR data (Panel g) of cell cycle sort experiments were analyzed. Error bars show standard deviation. Asterisks mark statistical significance (p < 0.05)

Full size image

Gene expression changes observed by qRT-PCR experiments of six genes chosen upon microarray analysis confirmed the microarray results in all of the three cells (Fig. 2, Panel c-e). Specifically, all the six genes chosen for the qPCR validation were present in the significant (NCI-H295R, HeLa) or FC > 2 (HDFa) lists. Moreover, ARHGAP11A, KIF14 and GTSE1 were previously found to be expressed in a cell cycle dependent fashion in primary fibroblasts [5] and HeLa cells [4], while ASPM and SKA1 genes were found to be cell cycle regulated in primary fibroblasts [5]. The successful validation of these well-known cell cycle genes in all three cell types analyzed here further confirms our cell cycle sorting method.

Functional bioinformatics analysis was used to detect altered pathways based on our microarray results. As a further confirmation of our method, ”Cell cycle” ”Cellular assembly and organization” and ”DNA replication, recombination and repair” were the molecular and cellular functions most concerned by gene expression changes in all three cells (Fig. 2, Panel f-h).

Comparison of cell cycle dependent expression between cell cycle sort and former synchronization based data

Several conflicting arguments arose on the applicability of synchronization procedures to define transcripts with cycling expression in unperturbed cells [7]. Therefore we aimed to compare expression changes between cell cycle phases detected by gene expression profiling in synchronization and cell cycle sort based experiments. Because synchronization based time course gene expression data in adrenocortical cell line have not been previously published, comparisons were made with primary fibroblasts and HeLa cells. Pearson’s method showed significant correlation between gene expression changes observed in synchronization based and cell cycle sort based experiments, confirming previous synchronization experiments by a synchronization-free method in unperturbed cells (Fig. 3, Panel a-c, Additional file 2: Figures S3 and S4, Additional file 1: Table S4).

Additionally, Gene Ontology (GO) Term analysis was performed on the HeLa cell cycle dependent transcriptional program to analyze the possible difference in biological processes affected by cell cycle sort and synchronization procedures. As both of cell cycle sort-based and synchronization-based results are only applicable in HeLa cells, we performed the analysis on three gene lists: genes unique to the HeLa cell cycle sort experiment (unique HeLa SORT), genes unique to the HeLa synchronization experiment (unique HeLa synchr) and the overlap between these two lists. All three lists were enriched with cell cycle-related processes; however, the overlap between the two experiments presented the most significant enrichment of cell cycle-associated biological processes, cross-validating important cell cycle genes detected by both the synchronization-based and cell cycle sort-based procedures. All the GO terms detected in the unique HeLa SORT list were detected in the overlap list, however, interestingly, five out of eight GO terms detected in the unique HeLa synchr list were unique to this list of genes, not being present in the analysis of the unique HeLa SORT or overlap gene lists (Table 1 and Additional file 1: Table S5).

Table 1 GO term analysis of the cell cycle dependent transcriptional program of HeLa cells
Full size table

Magnitude of gene expression alterations during cell cycle progression in untransformed and cancer cells

QRT-PCR validation of microarray experiments (Fig. 2, Panel c-e) indicated that gene expression changes might be characterized by different amplitudes in primary vs. cancer cells. Therefore, we analyzed the expression profiles and cell cycle dynamics of genes displaying altered expression between cell cycle phases in both primary untransformed (HDFa) and transformed cancer (HeLa) cells (127 genes present in both HDFa SORT and HeLa SORT gene lists, Fig. 3, Panel a). Significantly lower expression values were found in primary untransformed compared to cancer cells in G1, S and G2 phases as well (Fig. 3, Panel d-e). For the analysis of mRNA dynamics during cell cycle in untransformed and cancer cells, differences in mean fold changes of expression of genes commonly altered in HDFa and HeLa cell cycle sort experiments were calculated and evaluated (Fig. 3, Panel f-g). Among several significant alterations, a robust difference in mean fold change of gene expression was observed in G1/S transition between primary fibroblasts and cancer (NCI-H295R and HeLa) cells based on both microarray and qRT-PCR results. During the cell cycle, cycling genes had lower basal expression, but they demonstrated expression changes of significantly greater amplitude in primary non-transformed fibroblasts, than in transformed cancer (NCI-H295R and HeLa) cells.

High-throughput screening of miRNAs differently expressed during the human cell cycle

Three high-throughput platforms (microarray, TaqMan Low Density Array and Illumina small RNA Sequencing) of miRNA expression were used to detect cell cycle dependent miRNA expression (Fig. 4, Panel a-d, Additional file 2: Figure S5). Among them, microarray (Fig. 4, Panel a and c) displayed the lowest dynamic range and was unable to detect miRNAs of altered expression between cell cycle phases in HDFa and NCI-H295R cells. TaqMan Low Density array (Fig. 4, Panel b) performed on RNA isolated from sorted NCI-H295R cells detected 8 miRNAs of altered expression between cell cycle phases (among which only hsa-miR-10b, hsa-miR-128a and hsa-miR-890 had fold change values exceeding 2), however qRT-PCR validation of selected miRNAs failed to confirm the results (Additional file 2: Figure S6). Among the three platforms used in our study, small RNA sequencing was found to have the largest dynamic range in detection of miRNA expression alterations (Fig. 4, Panel d and Additional file 2: Figure S5). Still, statistical analysis detected only 11 miRNAs with altered expression in HeLa cells, of which only four miRNAs (hsa-miR-146b, hsa-miR-577, hsa-miR-877 and hsa-miR-193b*) had FC > 2 expression change between cell cycle phases. QRT-PCR measurements, similar to that of TLDA validation attempt, failed to validate differential expression in NCI-H295R and HeLa cells (Additional file 2: Figures S5 and S6). For further validation four other miRNAs showing stable expression in different cell types and cell cycle phases based on the high-throughput data were selected for qRT-PCR control analysis, which confirmed the stable expression pattern (Additional file 2: Figure S6).

Fig. 4
figure 4

Analysis of cell cycle dependent miRNA expression by high-throughput screenings methods and individual qRT-PCR in HDFa, NCI-H295R and HeLa cells. Panels a-d: Fold change (log2) of miRNA expression in S/G1, G2/S and G2/G1 phases observed by microarray (Panel a: NCI-H295R cells, Panel c: HDFa cells), TaqMan Low Density Array (Panel b: NCI-H295R cells) and Illumina Small RNA Sequencing (Panel d: HeLa cells). Grey background corresponds to fold change ≤ 2 between cell cycle phases. Colored squares represent miRNAs with significantly different expression in cell cycle phases obtained by statistical analysis of high-throughput data. Colored circles represent some members of the hsa-miR-16 family, if available. For further high-throughput data of NCI-H295R cells (Illumina Small RNA Sequencing) see Additional 2: Figure S5. For qRT-PCR measurement data see Additional 2: Figure S6. Panels e-g: QRT-PCR measurements of hsa-miR-16 family members hsa-miR-16, hsa-miR-15a and hsa-miR-503 in HDFa (Panel e), NCI-H295R (Panel f) and HeLa (Panel g) cells. In each case ΔCt values were normalized to G1 phase (ΔΔCt) and were subjected to FC = 2-ΔΔCt transformation. Error bars show standard deviation. Asterisks mark statistical significance (p < 0.05)

Full size image

Among several cell cycle regulator miRNAs, members of the hsa-miR-16 family were found to display dynamic changes in expression between serum-starved G0 and actively proliferating state [23]. Therefore, we analyzed expression changes of the hsa-miR-16 family members: hsa-miR-16, hsa-miR-15a and hsa-miR-503 in our high-throughput data (Fig. 4, Panel a-d and Additional file 2: Figure S5, Panel a) and performed qRT-PCR analysis as well (Fig. 4, Panel e-g). In the case of hsa-miR-15a, small RNA sequencing detected two-fold alteration in expression in NCI-H295R and HeLa cells as well and qRT-PCR analysis further confirmed some small, but significant expression changes in all three cell types.

Discussion

With the advent of high-throughput transcriptional profiling, the precise analysis of the cell cycle transcription program became possible. Time course gene expression data after various synchronization procedures were used for gene expression profiling on microarray platforms [4–6] which demonstrated the cell cycle transcription program in various cell types including untransformed primary [5, 6] and transformed cancer cells [4].

However, questions were raised on the applicability of synchronization procedures in the analysis of cell cycle dependent transcriptional profiling [7]. In fact, growth imbalance and unscheduled expression of key cell cycle factor cyclins were demonstrated due to synchronization procedures [8, 9]. Additionally, cells failed to retain their synchronization relatively soon after depriving the synchronization agent from growth medium [5, 10]. A rigorous set of criteria was introduced by Shedden and Cooper concerning the analysis of cell cycle dependent gene expression [7] and, therefore, novel methods satisfying these criteria were needed. Centrifugal elutriation [25] and cell cycle sort [26–28] are two synchronization-free methods satisfying the Shedden and Cooper criteria without perturbing the cell cycle machinery. However, to our knowledge no detailed study using these methods to determine cell cycle dependent gene expression in human cells has been published to date. Upon centrifugal elutriation, cells are differentiated upon the size of each cell [25], while during cell cycle sorting by FACS the amount of DNA in each cell is used for separation [27, 28]. Both methods achieved efficient separation of cell cycle phases [12, 13, 27–29]. Although proper segregation of cell cycle phases by cell cycle sort may be optimal in cells lacking aneuploidy, cells presenting aneuploidy may as well be introduced to cell cycle sort until proper further confirmation of cell cycle sort (detection of the different expression of key cell cycle regulators e.g. Tyr-15 phosphorylation of CDC-2 [24]) supports the results. In our results, dynamic changes of Tyr-15 phosphorylation of CDC-2 in HeLa show as great variability as in NCI-H295R cells, confirming successful segregation of cell cycle phases in HeLa cells as well.

Here, using the cell cycle sort method, we report a comparative analysis of the cell cycle dependent gene transcriptional profile of human untransformed primary (fibroblasts) and cancer (adrenocortical and cervical) cells. Functional bioinformatics analysis revealed cell cycle related molecular and cellular functions to be mostly concerned with these transcriptional alterations. As time-course gene expression data from previous synchronization-based studies regarding primary human fibroblasts [5] and HeLa cells [4] were available, we performed a re-analysis of these data to compare synchronization-based and cell cycle sorting-based results. Upon the significant correlation between our data and these earlier results, we conclude that data obtained from cell cycle sort experiments confirm earlier results demonstrating cell cycle dependent gene expression in human cells, as well as it satisfies the rigorous criteria described above [7].

Moreover, GO term analysis was performed to assess biological processes related to the HeLa cell cycle dependent transcriptional program. The overlap between HeLa cell cycle sort and HeLa synchronization experiments showed a robust enrichment of cell cycle-related GO terms, cross-validating the key players of the cell cycle dependent transcriptional program. Additionally, cell cycle-related biological processes were enriched in both the unique to HeLa cell cycle SORT and the unique to HeLa synchonization gene lists. However, interestingly, the majority of GO terms detected in the unique HeLa synchronization gene list were absent from the overlap gene list, indicating some specific mechanisms related to synchronization procedures. The specific presence of ”response to DNA damage stimulus” and ”cellular response to stress” and the induction of ”DNA repair” GO terms confirms the replication stress as a consequence of synchronization procedures [8, 9, 30]. These analyses further confirm that our synhronization-free method of cell cycle sort charactherises more specifically the cell cycle of unperturbed cells than the synchronization-based methods.

We also investigated the eventual difference of cell cycle regulated transcriptional program in human untransformed and cancer cells. Whitfield et al. demonstrated that genes exhibiting cell cycle regulated expression were overexpressed in malignant tumors reflecting the malignancy signature of neoplasms [4]. This was explained by the fact that tumors contain more cycling cells [4]. Cell cycle dynamics alter disproportionally during malignant transformation [31]: activation of oncogenes HRAS, SRC, MYC, CCND1, CCNE [32–34], and loss of tumor suppressor genes as PTEN [35] shortens G1 phase [31], while loss of key M phase regulators LZTS1 and LATS2 results in M phase shortening [31, 36, 37]. These alterations lead to a relatively larger portion of cells residing in S and G2 phases. Additionally, certain gene clusters were confirmed to exhibit cell cycle dependent expression in either primary untransformed or transformed cancer cells [5], differentiating cells upon malignant transformation.

Our results contribute to the notion of different transcriptional regulation in untransformed and cancer cells. Since we have analyzed only three human cell types of different tissue origin, we can not draw a definitive conclusion universal to the cell cycle effects of malignant transformation. However, based on our analysis we may hypothesize that genes displaying universal cell cycle dependent expression in untransformed and cancer cells display altered expression in each phase and dynamic changes of different amplitude (Fig. 5). MRNA expression was found to be higher in G1, S and G2 phases as well, therefore, in addition to altered cell cycle distribution, basal, phase-independent up-regulation of these cell cycle genes may as well constitute to the well observed higher expression in malignant cancers. Dynamic mRNA expression differences between G1 and S phases were of greater amplitude in untransformed primary cells than in cells undergoing malignant transformation. This may be explained by the longer and more tightly controlled G1 phase and G1/S transition observed in untransformed, primary cells [31], as it reflects the more precisely regulated cell cycle machinery in untransformed cells. Moreover, MYC amplification stimulates E2F expression in cancer cells, facilitating the commitment to cell division [38, 39]. This facilitated regulation of the G1/S transition may as well contribute to smaller expression changes of the cell cycle dependent transcriptional program.

Fig. 5
figure 5

Schematic presentation of the hypothesis concerning expression dynamics of mRNAs and miRNAs during the cell cycle phases in primary untransformed and cancer cells. Relative expression changes in G1, S and G2 phases of a representative cell cycle gene and miRNA in untransformed primary (Panel a) and transformed cancer (Panel b) cells. Relative expression changes of mRNAs exhibiting cell cycle dependent expression (dark blue) and cell cycle associated miRNAs (green) are shown. Blue double arrows mark the amplitude of dynamic expression range in cell cycle phases. Note the different phase-specific expression values and the different dynamic expression range in the two cell types. Upon our results, genes exhibiting cell cycle dependent expression profile in untransformed cells are characterized with lower expression levels throughout the cell cycle, however they possess greater variance in expression levels than in transformed, cancer cells. MiRNAs, however, do not display cell cycle dependent expression

Full size image

MiRNAs have a well established role in the regulation of the cell cycle [19]. Oncogenic (onco-miRs) and tumor-suppressor (TS-miRs) miRNAs were confirmed as modifiers of key cell cycle agents, accelerating or decelerating cell cycle progression [40, 41]. Long-term miRNA-mediated cell cycle changes contribute to malignant transformation in a variety of neoplasms [19, 42–46]. Additionally, the role of miRNA-mediated regulation has been confirmed in the transition from quiescent state to actively proliferating state [19, 23]. In particular, mitogenic stimuli enhances cell cycle progression by stimulating key transcriptional factors of the E2F family, which in turn enhances members of the hsa-let-7 and hsa-miR-16 families [17, 22, 23]. These well-known TS-miRs target key cell cycle cyclins as cyclin E, fine-tuning the proper cell cycle progression [17]. However, the proposed cell cycle dependent miRNA expression pattern [19] has not been thoroughly investigated, as according to our knowledge only one synchronization study identifying some cell cycle regulated miRNAs has been published to date [47].

Our study aiming to detect miRNA expression changes between cell cycle phases included the application of miRNA microarray, qPCR-based TLDA and Illumina small RNA sequencing. Microarrays are widely used for high throughput miRNA profiling and produce results which can be validated in high percentage by qPCR [48]. However, obtaining negative results prompted us for further analysis using qPCR-based TLDA and Illumina small RNA sequencing. The latter approach has larger dynamic range of detection allowing us to successfully detect smaller, but significant alterations [49, 50]. QPCR-based TLDA results however can be most successfully validated by single tube individual miRNA-specific qPCR as primer sequences used in TLDA does not differ.

Analysis of the three human cells (two cancer cell lines and one primary cell) on three high-throughput miRNA expression platform in our study revealed that miRNA expression profile throughout the cell cycle phases was quite stable (Fig. 5). Surprisingly, our systematic study using multiple high-throughput platforms indicated the lack of validable cell cycle dependent miRNA expression, and also showed that fold change differences are of small amplitude, especially in the light of the robust and explicit changes observed in mRNA expressions of the very same cell stage samples. More than 50 % of miRNA genes are located in cancer-associated genomic regions or in fragile sites, being continuously downregulated or deleted in cancers [40, 51]. Therefore, the loss of TS-miRs and the activation of onco-miRs are specifically involved in the long-term malignant transformation. Moreover, with the loss of genetic regions containing miRNA genes, the possibility of their dynamic regulatory functions throughout the cell cycle is lost as well [40, 51]. Additionally, miRNA-dependent gene regulation was found to be a much slower process than previously thought, due to certain bottlenecks related to the complex biogenesis and maturation processes or delays of miRNAs loading into Argonaute proteins [52]. Accelerated turnover was proposed to be necessary for certain miRNAs to be possibly involved in dynamic cell cycle regulation [52]. Such accelerated turnover has been confirmed in the case of the hsa-miR-16 family [23]. This family has been identified as a cluster of TS-miRs [19], downregulated in various types of cancers [19, 42–44] and associated with quiescent state [23]. Our results indicated some minor miRNA expression changes, especially in the case of hsa-miR-15a, however, expression changes were not fully congruent in the three cell types studied, suggesting that the well known, cell type specific expression of miRNAs may contribute to this phenomenon.

Finally, it is of utmost importance to address the limitations of our study. Firstly, upon the analysis of three human cell types of different tissue origin we can not draw a general conclusion concerning the difference of the expression dynamics of the universal cell cycle genes. Secondly, although different culture conditions may have some effect on the observed cell cycle differences, our method has much less effects on cell cycle alterations compared to previous serum shock-based or inhibitive synchronization-based processes.

Conclusions

In conclusion, successful utilization of cell cycle sort as a novel method for the analysis of cell cycle transcriptional program in our study confirmed the previously identified cell cycle transcriptional regulation. Different phase-dependent and phase-independent mRNA expression dynamics of cell cycle genes in human untransformed and cancer cells were revealed, reflecting the altered cell cycle machinery in cancer cells at the transcriptional level. Perhaps more interestingly, the application of various high-throughput platforms (microarray, TLDA, Small RNA Sequencing) for miRNA profiling showed that miRNA expression dynamics are unaltered during the active cell cycle at the G1/S and S/G2 transitions.

Methods

Cell culture

Human adrenocortical cancer cell line NCI-H295R and human cervical cancer cell line HeLa were obtained from the American Type Culture Collection (ATCC), while human dermal fibroblast (HDFa) cells were obtained from Gibco (Life Technologies). NCI-H295R cells were cultured in Dulbecco’s modified Eagle’s medium/Nutrient Mixture F-12 Ham (DMEM: F12) supplemented with 6.25 ng/ml insulin, 6.25 ng/ml transferrin, 6.25 ng/ml sodium selenite, 1.25 mg/ml bovine serum albumine, 5.35 ng/ml linoleic acid, 1 % HEPES, 1 % Penicillin-Streptomycin, 2.5 % L-glutamine (Sigma-Aldrich Chemical Co.) and 2.5 % Nu-Serum (BD Biosciences). HeLa cells were cultured in Dulbecco’s modified Eagle’s medium/Nutrient Mixture F-12 Ham (DMEM: F12, Sigma-Aldrich Chemical Co.) supplemented with 10 % fetal bovine serum (Gibco by Life Technologies) and 1 % antibiotic-antimycotic solution (Sigma-Aldrich Chemical Co.). HDFa cells were cultured in Medium 106 supplemented with low serum growth supplement (LSGS, Gibco by Life Technologies). All cells were cultured at 37 °C in a humidified 5 % CO2 atmosphere.

DNA content based fluorescence activated cell sorting (FACS)

HDFa, NCI-H295R and HeLa cells were cultured in 150 cm2 cell culture flasks until 90 % confluency. Cells were trypsinized, washed, resuspended in complete medium and counted. Vybrant DyeCycle Orange (Molecular Probes by Life Technologies) was used to stain genomic DNA stoichiometrically in living cells (approximate fluorescence excitation and emission maxima were 519 nm and 563 nm, respectively), and was added in 1:500 dilution to 1 × 106 cells/ml cell suspension. After incubation at 37 °C for 30 min, protected from light, cells were centrifuged at 1000 rpm for 10 min and were resuspended in the sort medium (Hank’s Balanced Salt Solution without Ca2+ and Mg2+, containing 2 % fetal calf serum). FACSAria III cell sorter (Becton-Dickinson, Franklin Lakes, NJ, USA) was used for cell cycle analysis and sorting using 488 nm Argon laser. The fluorescence emission of Vybrant DyeCycle Orange was separated by a 556 longpass filter and detected through a 585/42 bandpass filter. At least 100,000 events were collected for analysis. Upon cell cycle analysis, cell populations resembling G1, S and G2 phases were gated according to cellular DNA quantity. Sorting did not exceed 30 min and all sorted populations were validated by flow cytometry analysis. Data were analyzed by BD FACSDiva v6.1.3 software (BD Biosciences, San Jose, CA, USA). Thereafter, cells were centrifuged, washed with ice-cold PBS and resuspended in QIAzol lysis reagent (Qiagen) or Western blot lysis buffer for subsequent RNA or protein isolation, respectively. Until RNA isolation or Western blot, samples were stored at −80 °C. Optimization from the protocol supplied by the manufacturer included the use of sort medium, concentration of the cell suspension before FACS analysis, respecting an upper time limit for sorting and immediate FACS reanalysis upon every cell cycle-sorted population.

Protein isolation and Western blot

Samples were thawed on ice, sonicated with ultrasound and incubated on ice for 30 min. Thereafter, samples were centrifuged at 13000 rpm and 2 °C for 15 min. Protein concentration was determined according by Bradford method using Varioskan Flash spectral scanning reader (Thermo Scientific) [53]. Optical density was determined at 595 nm. Samples were mixed with β-mercaptoethanol containing Laemli buffer and were incubated at 99 °C for 5 min. Thereafter, equal amount of samples were loaded on a 10 % polyacrylamide gel and electrophoresis was conducted on a Mini Protean electrophoresis equipment (Bio-Rad). Overnight blotting at 4 °C was performed to transfer proteins to a PVDF membrane (Millipore, Billerica, MA). Blotting efficiency was determined by Ponceau staining. Membranes were blocked with 5 % non-fat dry milk in TBS for 60 min at room temperature, and were incubated with primary phospho-CDC-2 (Tyr15) antibody (Cell Signaling Technology, cat. No.: 9111, dilution: 1:500) at 4 °C for 16 h. Thereafter, membranes were washed 5 times with 0.05 % Tween-20 containing TBS, and were incubated with secondary antibody (Cell Signaling Technology, cat. No.: 7074, dilution: 1:2000). All antibodies were diluted in 1 % non-fat dry milk containing TBS. After exposure to SuperSignal West Pico Chemiluminescent Substrate (Thermo Scientific), signals were visualized by Kodak Image Station 4000MM Digital Imaging System. Thereafter membranes were stripped with mild stripping buffer (0.2 M glycine, 0.1 % sodium dodecyl sulfate, 0.1 % Tween-20, pH = 2.2) by gentle agitation for 45 min at room temperature, and were blocked again for subsequent detection of loading control β-actin (Cell Signaling Technology, cat. No.: 4967, dilution: 1:2000). Membrane blocking, antibody incubations and signal detection were carried out exactly as in the case of phospho-CDC-2 detection. Densitometry of the detected bands was performed by Kodak Image Station. β-actin was used as loading control.

RNA isolation, messenger RNA profiling and validation by qRT-PCR

Total RNA was isolated using miRNeasy Mini Kit (Qiagen), according to the manufacturer’s instructions and was eluted in 50 uL nuclease-free water (Qiagen). RNA concentration and integrity was determined by the Agilent Bioanalyzer 2100 system (Agilent Technologies, Additional file 2: Figure S1, Additional file 1: Table S1).

High throughput profiling of gene and miRNA expression

Gene expression profiling

Gene expression profiling was performed on 100 ng RNA isolated from sorted G1, S and G2 phases of HDFa, NCI-H295R and HeLa cells. In all, 24 samples (2 or 3 samples of each phase) were analyzed using Agilent whole human genome 4x44K microarray slides (Agilent Technologies) following the manufacturer’s protocol [54].

miRNA expression profiling using microarray

MiRNA expression profiling was performed on 100 ng RNA isolated from sorted G1, S and G2 phases of HDFa and NCI-H295R cells. In all, 16 samples (2 or 3 samples of each phase) were analyzed. The miRNA expression profiling using microarray followed the manufacturer’s protocol [48]. Total RNA was labeled with Cy3 and amplified using Low Input Quick Amp Labeling Kit according to the manufacturer’s instructions. After RNA purification, labeled RNA was hybridized to Agilent 8 × 15 K Human miRNA Microarray Release 12.0. slides (Agilent Technologies), according to the manufacturer’s instructions. After washing, array scanning and feature extraction was performed by Agilent DNA Microarray Scanner and Feature Extraction Software 11.0.1.

miRNA expression profiling using TaqMan Low Density Array (TLDA)

RNA isolated from two samples of G1 and three samples of S and G2 phase-sorted NCI-H295R cells were studied using TaqMan Low Density Array (TLDA) cards, according to the manufacturer’s instructions. The miRNA expression profiling using TLDA was performed as previously reported [55]. 30 ng of total RNA was reverse transcribed and pre-amplified using Megaplex RT primer pool A and B and Megaplex PreAmp primers, respectively. Quantitative real-time PCR were carried out in TaqMan Human MicroRNA Array A and B on a 7900HT Real time PCR System (Applied Biosystems by Life Technologies).

miRNA expression profiling using Illumina small RNA sequencing

Two samples of each cell cycle phase of HeLa cells (six samples) and one sample of each cell cycle phase of pooled sorted NCI-H295R cells (threee samples) were analyzed. Small RNA sequencing was performed at BGI using Illumina Small RNA Sequencing Platform. For library preparation TruSeq Small RNA library preparation kit (Illumina, San Diego, California) was used. Sequencing was performed by SE50 with Illumina HiSeq2000, and 10 Mb clean reads were analyzed followed by routine algorithms (BGI Tech Solutions, Tai Po, Hong Kong).

qRT-PCR validation

For the gene expression qRT-PCR experiments, 30 ng of total RNA was reverse transcibed using SuperScript VILO cDNA synthesis kit according to the manufacturer’s instructions (Applied Biosystems by Life Technologies). Gene expression was quantified using predesigned Taqman probes (Additional file 1: Table S2, Applied Biosystems by Life Technologies) on a 7500 Fast Real-time PCR system (Applied Biosystems by Life Technologies). Gene expression data were normalized to the relative expression of ACTB.

For the miRNA expression qRT-PCR experiments, 5 ng of total RNA was reverse transcribed and quantified using TaqMan microRNA reverse transcription kit (Applied Biosystems by Life Technologies) and predesigned TaqMan probes (Additional file 1: Table S2, Applied Biosystems by Life Technologies) on a 7500 Fast Real-time PCR system (Applied Biosystems by Life Technologies). MiRNA expression data were normalized to the relative expression of RNU48.

All measurements were performed in triplicate (three biological, two technical replicates). Expression level was calculated by the ΔCt(S-phase) – ΔCt(G1-phase) and the ΔCt(G2-phase) – ΔCt(G1-phase) (ΔΔCt) methods.

Pathway analysis

Ingenuity Pathway Analysis (IPA, Ingenuity Systems) was used to detect molecular and cellular functions altered between cell cycle phases. Δ(G2-G1) gene expression changes of significantly differently expressed genes (NCI-H295R and HeLa) or genes with fold change > 2 expression (HDFa) were subjected to IPA core analysis.

Microarray data from former studies

Two former microarray studies identifying cell cycle dependent expression of mRNA transcripts in human primary fibroblasts [5] and HeLa cells [4] using synchronization based procedures were selected. Processed data from these experiments were downloaded from http://genome-www.stanford.edu/Human-CellCycle/HeLa/ [4] and from the European Bioinformatics Institute Array Express database (http://www.ebi.ac.uk/arrayexpress/experiments/E-TABM-263/) [5] and were re-analyzed [4, 5]. Upon published FACS analysis data, time points with highest levels of synchronous populations of each cell cycle phase were chosen to represent G1, S and G2 phases, respectively. Difference in gene expression between phases was calculated upon difference of normalized expression of a certain gene between time points representing each phase.

In these comparisons only those gene expression alterations were used where the cell cycle sort indicated cell cycle dependent gene expression changes (FC > 2 genes of HDFa and significant genes of HeLa experiment).

Gene ontology term analysis

Gene Ontology (GO) term analysis was performed to detect biological processes with enriched genes in the HeLa cell cycle transcriptional program. The online functional annotation tool of DAVID Bioinformatics Resources version 6.7 (https://david.ncifcrf.gov/) with Gene Ontology for biological processes (category: GOTERM_BP_FAT) was used. The input gene lists for the analysis were the genes unique to HeLa SORT experiment (HeLa SORT \ HeLa synchr), unique to HeLa synchronization experiment (HeLa synchr \ HeLa SORT) and the overlap between these two lists (HeLa SORT ∩ HeLa synchr). Bonferroni-corrected p-values < 0.05 were considered statistically significant [4].

Analysis of gene expression levels and cell cycle dynamics in primary and cancer cells

Analysis of gene expression levels and cell cycle dynamics in different cell types were performed by investigating the changes of 127 genes found to be cell cycle dependently expressed in HDFa and HeLa cells. Upon combined normalization of all cell cycle sort-based gene expression microarrays, normalized intensity values in each cell type in each cell cycle phase were compared. For the analysis of the gene expression dynamism during cell cycle progression the absolute values of fold changes between cell cycle phases were calculated and were subjected to comparison between HDFa, NCI-H295R and HeLa cell types.

Results of qRT-PCR experiments in 10 (Additional file 1: Table S2) out of these 127 genes were also subjected to these analyses. ΔCt values normalized to ACTB expression and absolute values of fold changes in cell cycle phases were calculated and were compared in all cell types.

Statistical analysis

Statistical analysis of the microarray data was performed by GeneSpring 12.6 (Agilent Technologies) software. Total signal normalization at the 75th percentile of raw signal values and baseline transformation at the median of all samples following Agilent’s recommendation were performed. Differently expressed genes between G1, S and G2 phases were detected by one-way ANOVA followed by Tukey’s Honestly Significant Difference post hoc test and Benjamini-Hochberg correction for multiple measurements.

ΔCt levels of individually measured mRNA and miRNA transcripts obtained by qRT-PCR measurements and subsequent normalization to housekeeping transcripts (ACTB or RNU48) were subjected to Students’ two sided independent samples T-test. Differences were analyzed between G1-S, S-G2 and G1-G2 phases, respectively. Center values shown are the average of replicate experiments.

On genes displaying cell cycle dependent expression revealed by cell cycle sort, Pearson’s correlation was used to calculate correlation between expression changes detected by different (cell cycle sort and various synchronization) methods.

Student’s two-sided paired samples T-test was used to detect difference in normalized expression of genes expressed in a cell cycle dependent manner between various cell types. Student’s two-sided independent samples T-test was used to detect difference in absolute values of fold change of cell cycle dependently expressed genes of various cell types. In all comparisons p-value <0.05 was considered statistically significant.

Statistical analysis for miRNA expression analysis of TLDA card was performed using Real-Time StatMinerTM software (Integromics, Granada, Spain). Expression level was calculated by the ΔΔCt method, and fold changes were obtained using the formula 2-ΔΔCt. Following quality control, expression levels were normalized to the geometric mean of all expressed miRNAs. One-way ANOVA was used to detect significantly altered expression. In all comparisons p-value <0.05 was considered statistically significant.

For identification of differentially expressed miRNAs of Small RNA Sequencing experiments edgeR package version 3.8.6 in R was used. Alignment to MirBase version 21.0 mature miRNA database was performed on reads longer than 18 nucleotides with maximum 1 mismatch. The input data for edgeR package were the pair of phases (G1-S, S-G2, G1-G2) with two samples for each phases. The classical exact T-Test and TMM normalization were applied. Benjamini and Hochberg’s algorithm was used to control the false discovery rate (FDR). The difference was statistically significant when both the p-value and the FDR was <0.05.

References

  1. Cerqueira NM, Fernandes PA, Ramos MJ. Understanding ribonucleotide reductase inactivation by gemcitabine. Chemistry. 2007;13(30):8507–15. doi:10.1002/chem.200700260.

    Article  CAS  PubMed  Google Scholar 

  2. Nitiss JL. Targeting DNA, topoisomerase II in cancer chemotherapy. Nat Rev Cancer. 2009;9(5):338–50. doi:10.1038/nrc2607.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Schwartz GK, Shah MA. Targeting the cell cycle: a new approach to cancer therapy. J Clin Oncol Off J Am Soc Clin Oncol. 2005;23(36):9408–21. doi:10.1200/JCO.2005.01.5594.

    Article  CAS  Google Scholar 

  4. Whitfield ML, Sherlock G, Saldanha AJ, Murray JI, Ball CA, Alexander KE, et al. Identification of genes periodically expressed in the human cell cycle and their expression in tumors. Mol Biol Cell. 2002;13(6):1977–2000. doi:10.1091/mbc.02-02-0030.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Bar-Joseph Z, Siegfried Z, Brandeis M, Brors B, Lu Y, Eils R, et al. Genome-wide transcriptional analysis of the human cell cycle identifies genes differentially regulated in normal and cancer cells. Proc Natl Acad Sci U S A. 2008;105(3):955–60. doi:10.1073/pnas.0704723105.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Cho RJ, Huang M, Campbell MJ, Dong H, Steinmetz L, Sapinoso L, et al. Transcriptional regulation and function during the human cell cycle. Nat Genet. 2001;27(1):48–54. doi:10.1038/83751.

    Article  CAS  PubMed  Google Scholar 

  7. Shedden K, Cooper S. Analysis of cell-cycle-specific gene expression in human cells as determined by microarrays and double-thymidine block synchronization. Proc Natl Acad Sci U S A. 2002;99(7):4379–84. doi:10.1073/pnas.062569899.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Darzynkiewicz Z, Halicka HD, Zhao H, Podhorecka M. Cell synchronization by inhibitors of DNA replication induces replication stress and DNA damage response: analysis by flow cytometry. Methods Mol Biol. 2011;761:85–96. doi:10.1007/978-1-61779-182-6_6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Gong J, Traganos F, Darzynkiewicz Z. Growth imbalance and altered expression of cyclins B1, A, E, and D3 in MOLT-4 cells synchronized in the cell cycle by inhibitors of DNA replication. Cell Growth Differ Mole Biol J Am Assoc Cancer Res. 1995;6(11):1485–93.

    CAS  Google Scholar 

  10. Iyer VR, Eisen MB, Ross DT, Schuler G, Moore T, Lee JC, et al. The transcriptional program in the response of human fibroblasts to serum. Science. 1999;283(5398):83–7.

    Article  CAS  PubMed  Google Scholar 

  11. Tobey RA, Valdez JG, Crissman HA. Synchronization of human diploid fibroblasts at multiple stages of the cell cycle. Exp Cell Res. 1988;179(2):400–16.

    Article  CAS  PubMed  Google Scholar 

  12. Orlando DA, Lin CY, Bernard A, Wang JY, Socolar JE, Iversen ES, et al. Global control of cell-cycle transcription by coupled CDK and network oscillators. Nature. 2008;453(7197):944–7. doi:10.1038/nature06955.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Simmons Kovacs LA, Mayhew MB, Orlando DA, Jin Y, Li Q, Huang C, et al. Cyclin-dependent kinases are regulators and effectors of oscillations driven by a transcription factor network. Mol Cell. 2012;45(5):669–79. doi:10.1016/j.molcel.2011.12.033.

    Article  CAS  PubMed  Google Scholar 

  14. Valencia-Sanchez MA, Liu J, Hannon GJ, Parker R. Control of translation and mRNA degradation by miRNAs and siRNAs. Genes Dev. 2006;20(5):515–24. doi:10.1101/gad.1399806.

    Article  CAS  PubMed  Google Scholar 

  15. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–97.

    Article  CAS  PubMed  Google Scholar 

  16. Schultz J, Lorenz P, Gross G, Ibrahim S, Kunz M. MicroRNA let-7b targets important cell cycle molecules i7n malignant melanoma cells and interferes with anchorage-independent growth. Cell Res. 2008;18(5):549–57. doi:10.1038/cr.2008.45.

    Article  CAS  PubMed  Google Scholar 

  17. Ofir M, Hacohen D, Ginsberg D. MiR-15 and miR-16 are direct transcriptional targets of E2F1 that limit E2F-induced proliferation by targeting cyclin E. Mole Cancer Res MCR. 2011;9(4):440–7. doi:10.1158/1541-7786.MCR-10-0344.

    Article  CAS  Google Scholar 

  18. Butz H, Liko I, Czirjak S, Igaz P, Khan MM, Zivkovic V, et al. Down-regulation of Wee1 kinase by a specific subset of microRNA in human sporadic pituitary adenomas. J Clin Endocrinol Metab. 2010;95(10):E181–91. doi:10.1210/jc.2010-0581.

    Article  CAS  PubMed  Google Scholar 

  19. Bueno MJ, Malumbres M. MicroRNAs and the cell cycle. Biochim Biophys Acta. 2011;1812(5):592–601. doi:10.1016/j.bbadis.2011.02.002.

    Article  CAS  PubMed  Google Scholar 

  20. Johnson SM, Grosshans H, Shingara J, Byrom M, Jarvis R, Cheng A, et al. RAS is regulated by the let-7 microRNA family. Cell. 2005;120(5):635–47. doi:10.1016/j.cell.2005.01.014.

    Article  CAS  PubMed  Google Scholar 

  21. Dong Q, Meng P, Wang T, Qin W, Qin W, Wang F, et al. MicroRNA let-7a inhibits proliferation of human prostate cancer cells in vitro and in vivo by targeting E2F2 and CCND2. PLoS One. 2010;5(4):e10147. doi:10.1371/journal.pone.0010147.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Bueno MJ. Gomez de Cedron M, Laresgoiti U, Fernandez-Piqueras J, Zubiaga AM, Malumbres M. Multiple E2F7-induced microRNAs prevent replicative stress in response to mitogenic signaling. Mol Cell Biol. 2010;30(12):2983–95. doi:10.1128/MCB.01372-09.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Rissland OS, Hong SJ, Bartel DP. MicroRNA destabilization enables dynamic regulation of the miR-16 family in response to cell-cycle changes. Mol Cell. 2011;43(6):993–1004. doi:10.1016/j.molcel.2011.08.021.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Morla AO, Draetta G, Beach D, Wang JY. Reversible tyrosine phosphorylation of cdc2: dephosphorylation accompanies activation during entry into mitosis. Cell. 1989;58(1):193–203.

    Article  CAS  PubMed  Google Scholar 

  25. Leman AR, Bristow SL, Haase SB. Analyzing transcription dynamics during the budding yeast cell cycle. Methods Mol Biol. 2014;1170:295–312. doi:10.1007/978-1-4939-0888-2_14.

    Article  PubMed  Google Scholar 

  26. Henderson L, Bortone DS, Lim C, Zambon AC. Classic “broken cell” techniques and newer live cell methods for cell cycle assessment. Am J Physiol Cell Physiol. 2013;304(10):C927–38. doi:10.1152/ajpcell.00006.2013.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Juan G, Hernando E, Cordon-Cardo C. Separation of live cells in different phases of the cell cycle for gene expression analysis. Cytometry. 2002;49(4):170–5. doi:10.1002/cyto.10173.

    Article  PubMed  Google Scholar 

  28. Van der Aa N, Cheng J, Mateiu L, Zamani Esteki M, Kumar P, Dimitriadou E, et al. Genome-wide copy number profiling of single cells in S-phase reveals DNA-replication domains. Nucleic Acids Res. 2013;41(6):e66. doi:10.1093/nar/gks1352.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Archer SK, Inchaustegui D, Queiroz R, Clayton C. The cell cycle regulated transcriptome of Trypanosoma brucei. PLoS One. 2011;6(3):e18425. doi:10.1371/journal.pone.0018425.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Kurose A, Tanaka T, Huang X, Traganos F, Darzynkiewicz Z. Synchronization in the cell cycle by inhibitors of DNA replication induces histone H2AX phosphorylation: an indication of DNA damage. Cell Prolif. 2006;39(3):231–40. doi:10.1111/j.1365-2184.2006.00380.x.

    Article  CAS  PubMed  Google Scholar 

  31. Mizuno H, Nakanishi Y, Ishii N, Sarai A, Kitada K. A signature-based method for indexing cell cycle phase distribution from microarray profiles. BMC Genomics. 2009;10:137. doi:10.1186/1471-2164-10-137.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Liu JJ, Chao JR, Jiang MC, Ng SY, Yen JJ, Yang-Yen HF. Ras transformation results in an elevated level of cyclin D1 and acceleration of G1 progression in NIH 3 T3 cells. Mol Cell Biol. 1995;15(7):3654–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Wimmel A, Lucibello FC, Sewing A, Adolph S, Muller R. Inducible acceleration of G1 progression through tetracycline-regulated expression of human cyclin E. Oncogene. 1994;9(3):995–7.

    CAS  PubMed  Google Scholar 

  34. Karn J, Watson JV, Lowe AD, Green SM, Vedeckis W. Regulation of cell cycle duration by c-myc levels. Oncogene. 1989;4(6):773–87.

    CAS  PubMed  Google Scholar 

  35. Sun H, Lesche R, Li DM, Liliental J, Zhang H, Gao J, et al. PTEN modulates cell cycle progression and cell survival b7y regulating phosphatidylinositol 3,4,5,-trisphosphate and Akt/protein kinase B signaling pathway. Proc Natl Acad Sci U S A. 1999;96(11):6199–204.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Vecchione A, Croce CM, Baldassarre G. Fez1/Lzts1 a new mitotic regulator implicated in cancer development. Cell Div. 2007;2:24. doi:10.1186/1747-1028-2-24.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Yabuta N, Okada N, Ito A, Hosomi T, Nishihara S, Sasayama Y, et al. Lats2 is an essential mitotic regulator required for the coordination of cell division. J Biol Chem. 2007;282(26):19259–71. doi:10.1074/jbc.M608562200.

    Article  CAS  PubMed  Google Scholar 

  38. Leung JY, Ehmann GL, Giangrande PH, Nevins JR. A role for Myc in facilitating transcription activation by E2F1. Oncogene. 2008;27(30):4172–9. doi:10.1038/onc.2008.55.

    Article  CAS  PubMed  Google Scholar 

  39. Dong P, Maddali MV, Srimani JK, Thelot F, Nevins JR, Mathey-Prevot B, et al. Division of labour between Myc and G1 cyclins in cell cycle commitment and pace control. Nat Commun. 2014;5:4750. doi:10.1038/ncomms5750.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Zhang B, Pan X, Cobb GP, Anderson TA. microRNAs as oncogenes and tumor suppressors. Dev Biol. 2007;302(1):1–12. doi:10.1016/j.ydbio.2006.08.028.

    Article  CAS  PubMed  Google Scholar 

  41. Shenouda SK, Alahari SK. MicroRNA function in cancer: oncogene or a tumor suppressor? Cancer Metastasis Rev. 2009;28(3–4):369–78. doi:10.1007/s10555-009-9188-5.

    Article  CAS  PubMed  Google Scholar 

  42. Calin GA. D7umitru CD, Shimizu M, Bichi R, Zupo S, Noch E et al. Frequent deletions and down-regulation of micro- RNA genes miR15 and miR16 at 13q14 in chronic lymphocytic leukemia. Proc Natl Acad Sci U S A. 2002;99(24):15524–9. doi:10.1073/pnas.242606799.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Bonci D, Coppola V, Musumeci M, Addario A, Giuffrida R, Memeo L, et al. The miR-15a-miR-16-1 cluster controls prostate cancer by targeting multiple oncogenic activities. Nat Med. 2008;14(11):1271–7. doi:10.1038/nm.1880.

    Article  CAS  PubMed  Google Scholar 

  44. Takeshita F, Patrawala L, Osaki M, Takahashi RU, Yamamoto Y, Kosaka N, et al. Systemic delivery of synthetic microRNA-16 inhibits the growth of metastatic prostate tumors via downregulation of multiple cell-cycle genes. Mole Therapy J Am Soc Gene Therapy. 2010;18(1):181–7. doi:10.1038/mt.2009.207.

    Article  CAS  Google Scholar 

  45. Diaz R, Silva J, Garcia JM, Lorenzo Y, Garcia V, Pena C, et al. Deregulated expression of miR-106a predicts survival in human colon cancer patients. Genes Chromosomes Cancer. 2008;47(9):794–802. doi:10.1002/gcc.20580.

    Article  CAS  PubMed  Google Scholar 

  46. Guo X, Guo L, Ji J, Zhang J, Zhang J, Chen X, et al. miRNA-331-3p directly targets E2F1 and induces growth arrest in human gastric cancer. Biochem Biophys Res Commun. 2010;398(1):1–6. doi:10.1016/j.bbrc.2010.05.082.

    Article  CAS  PubMed  Google Scholar 

  47. Zhou JY, Ma WL, Liang S, Zeng Y, Shi R, Yu HL, et al. Analysis of microRNA expression profiles during the cell cycle in synchronized HeLa cells. BMB Rep. 2009;42(9):593–8.

    Article  CAS  PubMed  Google Scholar 

  48. Tombol Z, Eder K, Kovacs A, Szabo PM, Kulka J, Liko I, et al. MicroRNA expression profiling in benign (sporadic and hereditary) and recurring adrenal pheochromocytomas. Modern Pathol Off J United States Canadian Acad Pathol Inc. 2010;23(12):1583–95. doi:10.1038/modpathol.2010.164.

    Article  Google Scholar 

  49. Stokowy T, Eszlinger M, Swierniak M, Fujarewicz K, Jarzab B, Paschke R, et al. Analysis options for high-throughput sequencing in miRNA expression profiling. BMC Res Notes. 2014;7:144. doi:10.1186/1756-0500-7-144.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Tam S, De Borja R, Tsao MS, McPherson JD. Robust global microRNA expression profiling using next-generation sequencing technologies. Lab Investig J Tech Meth Pathol. 2014;94(3):350–8. doi:10.1038/labinvest.2013.157.

    Article  CAS  Google Scholar 

  51. Calin GA, Sevignani C, Dumitru CD, Hyslop T, Noch E, Yendamuri S, et al. Human microRNA genes are frequently located at fragile sites and genomic regions involved in cancers. Proc Natl Acad Sci U S A. 2004;101(9):2999–3004. doi:10.1073/pnas.0307323101.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Hausser J, Syed AP, Selevsek N, Van Nimwegen E, Jaskiewicz L, Aebersold R, et al. Timescales and bottlenecks in miRNA-dependent gene regulation. Mol Syst Biol. 2013;9:711. doi:10.1038/msb.2013.68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Bradford MM. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Anal Biochem. 1976;72:248–54.

    Article  CAS  PubMed  Google Scholar 

  54. Tombol Z, Szabo PM, Molnar V, Wiener Z, Tolgyesi G, Horanyi J, et al. Integrative molecular bioinformatics study of human adrenocortical tumors: microRNA, tissue-specific target prediction, and pathway analysis. Endocrine-Related Cancer. 2009;16(3):895–906. doi:10.1677/ERC-09-0096.

    Article  CAS  PubMed  Google Scholar 

  55. Butz H, Liko I, Czirjak S, Igaz P, Korbonits M, Racz K, et al. MicroRNA profile indicates downregulation of the TGFbeta pathway in sporadic non-functioning pituitary adenomas. Pituitary. 2011;14(2):112–24. doi:10.1007/s11102-010-0268-x.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

This work was supported by the Hungarian Academy of Sciences “Lendület” grant awarded to Attila Patocs (Lendület 2013), by the Hungarian Scientific Research Fund (OTKA, PD100648 (AP)) and by the Technology Innovation Fund, National Developmental Agency (KTIA-AIK-2012-12-1-0010).

Authors’ contributions

VKG participated in the design, laboratory work (FACS, expressional profiling and validation, Western blot), performed data analysis and interpretation and drafted the manuscript. EAT carried out FACS analysis and helped to draft the manuscript. KB performed Western blot analysis and helped to draft the manuscript. IL performed TLDA profiling and statistical analysis. OD performed statistical analysis on small RNA Seq data and helped to draft the manuscript. IK participated in the design and in Western blot analysis. JM participated in the design and FACS analysis and helped to draft the manuscript. KR participated in the design and helped to draft the manuscript. AP conceived of the study, participated in its design and coordination and helped to draft the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Author information

Authors and Affiliations

  1. 2nd Department of Medicine, Semmelweis University, Szentkirályi utca 46, 1088, Budapest, Hungary

    Vince Kornél Grolmusz & Károly Rácz

  2. “Lendület” Hereditary Endocrine Tumours Research Group, Hungarian Academy of Sciences, Semmelweis University, Szentkirályi utca 46, 1088, Budapest, Hungary

    Vince Kornél Grolmusz, István Likó, Ottó Darvasi & Attila Patócs

  3. Department of Immunology, Eötvös Loránd University, Pázmány Péter sétány 1/C, 1117, Budapest, Hungary

    Eszter Angéla Tóth & János Matkó

  4. 1st Department of Pathology and Experimental Cancer Research, Semmelweis University, Üllői út 26, 1085, Budapest, Hungary

    Kornélia Baghy & Ilona Kovalszky

  5. Molecular Medicine Research Group, Hungarian Academy of Sciences – Semmelweis University, Szentkirályi utca 46, 1088, Budapest, Hungary

    István Likó, Ottó Darvasi, Károly Rácz & Attila Patócs

  6. Department of Laboratory Medicine, Semmelweis University, Nagyvárad tér 4, 1089, Budapest, Hungary

    Attila Patócs

Authors
  1. Vince Kornél Grolmusz
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Eszter Angéla Tóth
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Kornélia Baghy
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. István Likó
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Ottó Darvasi
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Ilona Kovalszky
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. János Matkó
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Károly Rácz
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Attila Patócs
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Attila Patócs.

Additional files

Additional file 1: Table S1.

Characterization of cell cycle sorted cells and isolated RNA quantity in various cell types. Purity of cell cycle sort was determined by re-analyzing the sorted populations by FACS analysis: percentage was determined by the portion of cells residing in the gate previously designated for a certain cell cycle phase. Data of four replicate experiments. Data are given as × 1000 cells (sorted cells) and are shown as mean ± standard deviation. Table S2. Name and details of primers used for mRNA and miRNA expression qRT-PCR measurements. mRNA primers (Panel A, cat. No: 4331182) and miRNA primers (Panel B, cat. No: 4427975). All primers were from Applied Biosystems by Life Technologies. Table S3. Normalized expression of differently expressed genes between cell cycle phases in various cell types. Normalized expression of genes with fold change > 2 between cell cycle phases detected in HDFa cells (Panel A). Normalized expression of significantly differently expressed genes in cell cycle phases detected in NCI-H295R (Panel B) and HeLa (Panel C) cells. Note: For Panel B and C, genes are listed in the manner as shown in the heat map (Fig. 2, panel A and B, respectively). Table S4. List of genes shown on Venn diagram (Fig. 3, panel a). Genes are marked with “1” if being found cycling by either method (HDFa SORT, PF synchr, HeLa SORT, HeLa synchr). Gene IDs are Gene Symbols, if available or probe IDs. Table S5. List of HeLa cell cycle genes being present in enriched GO terms. Table 1 presents GO terms which are enriched in gene lists unique to HeLa SORT, HeLa synchronization experiments and the overlap beween HeLa SORT and synchronization lists. Gene symbols of genes being present in the gene lists and the enriched GO terms are shown here (HeLa SORT \ HeLa synchr - Panel A; HeLa synchr \ HeLa SORT – Panel B; HeLa SORT ∩ HeLa synchr – Panel C). (XLS 743 kb)

Additional file 2: Figure S1.

Quality and integrity of RNA isolated from cell cycle sorted populations of G1, S and G2 phases of HDFa, NCI-H295R and HeLa cells. Isolated RNA was analyzed on Agilent Bioanalyzer 2100 System. Representative results of G1, S and G2 populations sorted from HDFa (Panel A), NCI-H295R (Panel B) and HeLa (Panel C) cells are shown. RNA integrity number (RIN) was calculated if RNA concentration exceeded 10 ng/uL, therefore, in samples with lower RNA concentration RIN is not shown. Figure S2. Supplementary functional bioinformatics analysis of molecular and cellular functions concerned by gene expression alterations in cell cycle phases. Panel A: HDFa, Panel B: NCI-H295R, Panel C: HeLa. In Fig. 2, panels f-h only the five most significantly concerned networks are shown. Additional molecular and cellular functions are shown here. Figure S3. Pearson’s correlation analysis of gene expression changes in cell cycle phases of cell cycle sort and synchronization method in primary fibroblasts. Data of synchronization experiments: [5]. Synchronization methods: SS – serum starvation, SST – serum starvation followed by thymidine block. Correlation coefficients are shown. Asterisks mark statistical significance (p < 0.05). Figure S4. Pearson’s correlation analysis of gene expression changes in cell cycle phases of cell cycle sort and synchronization method in HeLa cells. Data of synchronization experiments: [4]. Synchronization methods: DT – double thymidine block, TN – thymidine followed by nocodazole block. Correlation coefficients are shown. Asterisks mark statistical significance (p < 0.05). Figure S5. Illumina Small RNA Sequencing and qRT-PCR measuremnts of cell cycle sorted NCI-H295R cells. Panel A: Fold change (log2) of miRNA expression in S/G1, G2/S and G2/G1 phases observed by Illumina Small RNA Sequencing. One pooled sample of each cell cycle phase was sequenced and, therefore, no statistical analysis has been performed. Grey background corresponds to fold change ≤ 2 between cell cycle phases. Colored circles represent some members of the hsa-miR-16 family, if available. Colored triangles represent two chosen miRNAs with FC > 2 expression between G2/G1 phases chosen to validate high-throughput data. Panel B: QRT-PCR measurements of 4 miRNAs chosen to validate Illumina Small RNA Sequencing data. Hsa-miR-202 and hsa-miR-132 were chosen upon FC > 2 expression in cell cycle phases. Hsa-miR-24-2* and hsa-miR-577 were chosen as G1 phase lacked miRNA expression, while in G2 phase miRNA expression was detectable. However, neither statistical examination detected significant difference (p < 0.05). Error bars show standard deviation. Figure S6. QRT-PCR measurements of miRNAs chosen to validate high-throughput data in HDFa, NCI-H295R and HeLa cells. Hsa-miR-10b, hsa-miR-128a and hsa-let-7 g were chosen upon TLDA screening of NCI-H295R cells. Hsa-let-7a, hsa-let-7e and hsa-let-7f were chosen upon Small RNA Sequencing of cell cycle sorted HeLa cells. MicroRNAs of stable expression (hsa-let-7i, hsa-miR-21, hsa-miR-22 and hsa-miR-222) in all three cell types and on all platforms were chosen to validate stable expression. Data from HDFa cells (Panel A), NCI-H295R cells (Panel B) and HeLa cells (Panel C). Additionally, hsa-miR-577 was also measured in HeLa cells, but it showed no cell cycle dependent expression pattern (Panel D). Error bars show standard deviation. Neither statistical examination detected significant difference (p < 0.05). (TIF 636 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Grolmusz, V.K., Tóth, E.A., Baghy, K. et al. Fluorescence activated cell sorting followed by small RNA sequencing reveals stable microRNA expression during cell cycle progression. BMC Genomics 17, 412 (2016). https://doi.org/10.1186/s12864-016-2747-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-016-2747-6

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us