Gene-expression profiles of microtubule perturbation in quiescent cells

To identify interphase-specific responses to microtubule damage, we first established conditions of healthy quiescence in a reference cell type for cell biology studies. Retinal pigment epithelial 1 (RPE1) cells, immortalized by the expression of human telomerase reverse transcriptase (hTert), were grown to and maintained as confluent cultures for 5 days prior to treatment with microtubule-destabilizing and stabilizing drugs. Cell-cycle profiling showed strong G1 arrest and no difference in cell-cycle state between control and drug treated cells after 24 h (S1A Fig). Combretastatin A-4 (CA4), which binds tubulin at the colchicine site, was chosen as a reference destabilizer, and paclitaxel (PTX) as a reference stabilizer [19] (S1B–S1E Fig). We chose drug concentrations based on tubulin partitioning between polymer and dimer and counting growing microtubule plus-tips labeled with end-binding protein 1 (EB1) [20]. Low doses (1 nM CA4, 3 nM PTX) were chosen to cause barely detectable effects, and high doses (100 nM CA4, 300 nM PTX) to cause complete loss of polymer (destabilizer) or complete loss of soluble dimer (stabilizer) (S1B–S1E Fig). We then compared gene-expression responses by RNA sequencing (RNA-seq) of polyadenylated (polyA+) mRNAs at 6 and 24 h post-drug treatment. These data revealed differential expression of multiple genes under all conditions relative to basal expression (Benjamini–Hochberg corrected p-value < 0.05, and >50 mRNA counts per million reads, Fig 1A–1D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Microtubule damage triggers DGE.

(A-D) DGE in RPE1 hTert cells treated with microtubule poisons for 6 (a and b) or 24 h (c and d). Data are represented as average Log₂FCs from three independent biological replicates, relative to basal expression in DMSO-treated control cells. Plotted are expression profiles of genes that have p-value < 0.05 and >50 mRNA reads per million, in 1 nM CA4- and 3 nM PTX-treated cells (a and c) and 100 nM CA4- and 300 nM PTX-treated cells (b and d). In red are depicted genes with FDR < 0.02 for low drug doses and < 0.001 for high drug doses, of which names are printed for top 20 genes. Plot quadrants are labeled from I–IV. (E and F) GSEA based on expression profiles of genes from quadrants II and IV, from cells treated with high drug doses for 6 (b) or 24 h (d). Plotted are −Log₁₀ BH-adjusted p-values of the top five enriched GO terms in the category Biological Process, with BH-adjusted p-value < 0.01. ACTB, β-actin; ACTN, actinin; ANKRD, Ankyrin Repeat Domain; BH, Benjamini–Hochberg; CA4, combretastatin A-4; CCDC, Coiled-coil domain containing; COX, Cytochrome C Oxidase; CTGF, Connective Tissue Growth Factor; CYR, Cytochrome Reductase; CYTB, Cytochrome B; DEG, differentially expressed gene; DGE, differential gene expression; EGR, Early Growth Response; ERRF, Endoplasmatic Reticulum-Related Factor; FDR, false discovery rate; GO, gene ontology; GSEA, gene set enrichment analysis; hTert, human telomerase reverse transcriptase; KHDRBS, KH RNA Binding Domain Containing, Signal Transduction Associated; LBH, Limb Bud And Heart Development; Log₂FC, Log₂ Fold Change; MALAT, Metastasis-Associated Lung Adenocarcinoma Transcript; MAP, microtubule-associated protein; MEG, Maternally Expressed Gene; ND, NADH Dehydrogenase; NEAT, Nuclear Paraspeckle Assembly Transcript; NKTR, Natural Killer Cell Triggering Receptor; OGT, O-Linked N-Acetylglucosamine Transferase; PTX, paclitaxel; RPE1, retinal pigment epithelial 1; SCG, Secretogranin; SLC, Solute Carrier Family; SRF, Serum Response Factor; SRM, Spermidine Synthase; THBS, Thrombospondin; TMX, Thioredoxin Related Transmembrane Protein; TM4SF, Transmembrane 4 L Six Family Member; TNFRSF, Tumor Necrosis Factor Receptor Superfamily Member; TTC, Tetratricopeptide Repeat Domain; TUBA, α-tubulin; TUBB, β-tubulin; XIST, X Inactive Specific Transcript.

https://doi.org/10.1371/journal.pbio.3000225.g001

DGE relative to vehicle treated cells (Fig 1A–1D) is presented to reveal genes that respond similarly to both drugs (quadrants I and III) versus differentially (quadrants II and IV, Fig 1A–1D). As expected, high drug doses caused extensive gene-expression changes. Surprisingly, subthreshold doses also caused many significant changes, showing that quiescent cells can detect even mild perturbation of microtubule dynamics. Based on high-drug-dose–expression profiles, we performed gene set enrichment analysis (GSEA), finding significantly enriched gene ontology (GO) terms (Fig 1E and 1F). The most prominently enriched are genes involved in G1/S transition of the mitotic cell cycle, which were up-regulated in CA4-treated cells and down-regulated in PTX-treated cells (Fig 1F). These data are consistent with previous reports that microtubule destabilization promotes DNA replication and cell-cycle re-entry, while stabilization inhibits it [21–23]. We further found highly enriched genes involved in apoptosis, wound healing, and cell migration, which were also up-regulated upon microtubule destabilization and down-regulated upon stabilization (Fig 1E and 1F). Surprisingly, we found only one cluster of microtubule-related genes involved in cytoskeleton-based transport (6 h and 24 h post-treatment) that were differentially regulated by microtubule stabilization versus destabilization (Fig 1E and 1F). We conclude that cells mount both differential and common responses to microtubule stabilization versus destabilization. Furthermore, they clearly detect damage that is subthreshold by biochemical tubulin-partitioning criteria. Notably, we did not observe changes in mitosis-related genes in quadrants I and III, consistent with the lack of mitotic arrest in our quiescent cultures. Some common DGE signatures were suggestive of general stress responses, such as apoptosis, and these deserve further analysis.

Tubulin genes stood out as the most differentially regulated genes in quadrant IV in all four plots, particularly in the low-dose perturbations (Fig 1A–1D). This observation was statistically confirmed by GSEA analysis, in which tubulin differential expression drove enrichment of GO terms such as “cytoskeleton-dependent intracellular transport” (Fig 1E and 1F). We conclude that coordinated change in multiple tubulin mRNAs was the strongest response to drug-induced microtubule damage, with destabilizing drugs decreasing tubulin mRNA concentrations and stabilizing drugs increasing them.

Differential regulation of most tubulin isoforms upon microtubule damage

To test if all tubulin mRNAs undergo similar regulation, we extracted the expression profiles of all detected tubulin genes in our data set, and their total mRNA counts (Fig 2A). We observed differential regulation of all highly expressed TUBA and TUBB mRNAs, but not centrosomal δ- (TUBD1) and ε-tubulins (TUBE1) (Fig 2A). Moreover, we found clear differential expression of centrosomal TUBG isoform 1 (TUBG1) and borderline regulation of isoform 2 (TUBG2) (Fig 2A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Microtubule damage triggers differential expression of all TUBA and TUBB isoforms.

(A, left panel) Expression profiles of all detected TUBA and TUBB isoforms in our DGE data set. Dendrogram on the left represents Pearson distance between expression profiles. Each column of the heatmap represents DGE in one treatment, labeled on the x-axis above and below the heatmap, relative to DMSO control. Each row represents a gene, labeled on the y-axis. In bold are differentially expressed genes that have p-value < 0.05. Color key is depicted in upper left corner. Data are represented as Log₂FC relative to DMSO control. (A, right panel) Average Log₂ mRNA CPMs in DMSO-treated control samples are presented as a bar chart. (B and C) Relative abundance of unspliced TUBA1A and TUBB pre-mRNA in control (DMSO) and cells treated with microtubule poisons (x-axis) for 6 (B) or 24 h (C). (D and E) Relative abundance of spliced TUBA1A and TUBB mRNA in control (DMSO) and cells treated with microtubule poisons (x-axis) for 6 (D) or 24 h (E). All the expression profiles are normalized to a reference gene (GAPDH or RPL19) and to DMSO control. Error bars in all panels represent standard deviation from three independent biological replicates. *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001 in paired Student t test compared to DMSO control. CA4, combretastatin A-4; CPM, count for each gene per million detected reads; DGE, differential gene expression; GAPDH, Glyceraldehyde 3-phosphate dehydrogenase; Log₂FC, Log₂ Fold Change; PTX, paclitaxel; RPL19, ribosomal protein L19; TUBA, α-tubulin; TUBB, β-tubulin; TUBD, δ-tubulin; TUBE, ε-tubulin; TUBG, γ-tubulin.

https://doi.org/10.1371/journal.pbio.3000225.g002

To generalize this finding, we reanalyzed two large, high-quality data sets deposited in the Gene Expression Omnibus (GEO) database that profiled DGE response to microtubule damage. In an extensive study that compared PTX with eribulin (ERB, a microtubule destabilizer) treatment of many breast, ovarian, and endometrial cancer cell lines [16], we confirmed differential regulation of all expressed TUBAs and TUBBs and TUBG1 (S2A Fig). Importantly, reanalyzing a study that compared the effect of microtubule destabilizers colchicine, vinblastine, and vincristine on rat heart endothelial cells [24], we show for the first time differential regulation of tubulin genes in vivo (GEO GSE19290, S2B Fig). We conclude that cells differentially regulate all the expressed TUBA and TUBB isoforms and TUBG1 upon microtubule damage, both ex vivo and in vivo.

The microtubule-damage–induced changes in tubulin mRNA concentrations that we observed were strongly suggestive of tubulin autoregulation, a post-translational gene-expression regulation mechanism [25]. RNA-seq of polyA+ mRNA does not distinguish between transcriptional and post-transcriptional regulatory mechanisms because the sample is enriched for spliced mRNA. Similarly, most microarray assays target exclusively the exonic sequences of mRNAs, making it impossible to distinguish the regulation of unspliced and spliced mRNA and draw conclusions about transcriptional versus post-transcriptional gene-expression regulation. To make this determination, we established a reverse-transcription quantitative PCR-based assay (RT-qPCR) to specifically measure transcriptional regulation through the expression levels of unspliced pre-mRNA and post-transcriptional regulation through the expression levels of spliced mRNA (S2C Fig). Using this approach, we measured two highly expressed tubulin genes, TUBA1A and TUBB, and two control housekeeping genes, Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) and ribosomal protein L19 (RPL19). We found no significant change in unspliced TUBA1A and TUBB pre-mRNA concentration in cells treated with CA4 or PTX (Fig 2B and 2C), showing that microtubule damage did not change the rate of tubulin gene transcription. However, levels of mature, spliced TUBA1A and TUBB mRNAs significantly diminished in CA4-treated cells and increased in PTX-treated cells (Fig 2D and 2E), consistent with our RNA-seq data. We conclude that post-transcriptional regulation of tubulin mRNA stability is the most prominent gene-expression response to microtubule damage. Importantly, we did not observe coregulation of any microtubule-interacting proteins, such as microtubule-associated, motor, or plus-tip–binding proteins. Thus, altered stability of microtubules only regulates the expression of tubulins, but not the other components of functional microtubules.

Bioinformatic analysis of the autoregulation signature reveals new microtubule biology

We next sought to investigate whether tubulin DGE is a general response to altered microtubule dynamics in conditions other than microtubule-targeted poisoning. The differential tubulin gene expression triggered by microtubule damage comprises a strong and specific signature that can be used to query publicly available DGE datasets in an unbiased manner and with the expectation of finding novel conditions that regulate microtubules. To test this approach, we used CLustering by Inferred Co-expression [26] (CLIC, https://gene-clic.org, Fig 3A)—a bioinformatic tool that mines approximately 3,500 publicly available human and mouse microarray studies deposited in the GEO database. Importantly, most of these studies are not designed to research cellular response to microtubule damage, providing an unbiased approach that can potentially reveal new microtubule biology.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Cells coordinate expression of TUBA and TUBB isoforms.

(A) Scheme of the bioinformatic approach. (B and C) Pearson expression correlation coefficients for a subset of abundantly expressed TUBA and TUBB isoforms across 417 human (B) and 122 mouse (C) publicly available Affymetrix chip data sets. In dark red are genes that show strong expression correlation (Pearson correlation coefficient = 1), and in gray are genes that show no expression correlation (Pearson correlation coefficient = 0). Color keys are represented to the left of each heatmap. High expression correlation for a subset of tubulin genes is marked as cluster CEM1. Tubulin genes that do not correlate in expression are marked as Null module. (D) Frequency diagram of groups of perturbations among the top 100 data sets with the highest tubulin gene-expression correlation from the human platform. CEM, Co-Expressed Modules; GEO, Gene Expression Omnibus; TUBA, α-tubulin; TUBB, β-tubulin; TUBD, δ-tubulin; TUBE, ε-tubulin; TUBG, γ-tubulin.

https://doi.org/10.1371/journal.pbio.3000225.g003

We first investigated whether tubulin genes were coordinately expressed across many sample and perturbation types. We found high Pearson expression correlation coefficients for the most broadly expressed tubulin genes in 417 human (Fig 3B) and 122 mouse studies (Fig 3C). More specialized and low-abundance tubulin isoforms displayed lower expression correlation coefficients. The obtained expression correlations confirm that cells coregulate the mRNAs for multiple TUBA and TUBB isoforms and TUBG1 in many cell and tissue types across a total of 539 different perturbations. These perturbations are candidates for novel conditions that regulate the microtubule cytoskeleton.

To understand what kinds of experimental perturbations induce coregulation of tubulin genes and potentially regulate the microtubule cytoskeleton, we rank-ordered the human data sets by descending Pearson correlation coefficient (S2 Table) and manually annotated the top 100 (Fig 3D and Table 1). Among high-ranked studies were multiple investigations of established microtubule drugs, such as PTX, or toxins known to perturb microtubules, such as bisphenol-A (Table 1). These served as positive controls that our bioinformatic analysis returned studies in which microtubules were perturbed. We also discovered many novel, to our knowledge, conditions, including virus infection, metabolite deprivation, exposure to industrial toxins, inflammation, and cell differentiation, which all tended to cause coordinated down-regulation of tubulin mRNAs (Fig 3D and Table 1). Surprisingly, only one process, oncogenic transformation, consistently up-regulated tubulin mRNAs (Fig 3D and Table 1). Hormone, vitamin, and chemical treatment, as well as perturbations of signaling, had molecule-specific effects on tubulin gene expression (Fig 3D, Table 1). In some studies, inspection of the CLIC report showed that coregulation of tubulin was not induced by the experimental perturbation in the title of the study, but rather by some other controlled and annotated experimental variable. For example, Steroid Receptor Coactivator-1 (SRC1) RNAi in A549 lung cancer cells did not change tubulin gene expression compared to control RNAi, but glucose withdrawal, which was included as a variable in the same study, suppressed it [27] (S3A Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. The top 100 data sets with highest tubulin-gene–expression Pearson correlation from the human platform.

Annotated are studies rank-ordered by decreasing tubulin-gene–expression Pearson correlation, with their associated GEO identifiers. Extracted from the CLIC report is differential tubulin gene expression. Unless NA, down-regulated (down) or up-regulated (up) tubulin mRNA levels are annotated. Based on published literature, perturbations in the listed studies are annotated into 11 groups, and their effect on microtubule cytoskeleton is marked as known (yes) or unknown (no).

https://doi.org/10.1371/journal.pbio.3000225.t001

To investigate whether our bioinformatic analysis returned new microtubule biology, we selected several perturbations, starting with nutrient deprivation, for validation by RT-qPCR and performed subsequent biochemical and microscopy-based inspection of the microtubule cytoskeleton. We confirmed that glucose and glutamine deprivation lower tubulin mRNA levels in cancer cell lines, further finding that this occurred by a mixture of transcriptional and mRNA stability mechanisms (see S3 Fig and related text). In both scenarios, nutrient deprivation caused destabilization of the microtubule cytoskeleton (S3D and 3H Fig). We conclude that various physiological and damaging inputs change tubulin gene expression either transcriptionally, thus affecting the microtubule cytoskeleton, or post-transcriptionally through changes in microtubule dynamics.

PI3K activity increases tubulin gene expression via microtubule stability and autoregulation

Multiple investigations of PI3K signaling scored highly in our bioinformatics search, including studies in which genetics were used to up-regulate the pathway activity and small-molecule inhibitors to down-regulate it. In each case, PI3K pathway activation correlated with an increase in tubulin mRNA levels (GEO GSE17785) and pathway inhibition with a decrease [28] (Table 1, S1 Table). One CLIC report suggested that MCF10A cells that expressed constitutively active PI3K with a mutation in the kinase domain [29] (H1047R) had increased tubulin mRNA levels (GEO GSE17785, Fig 4A). PI3K inhibition with GDC-0941 prevented this up-regulation of tubulin mRNA levels, suggesting involvement of PI3K in concerted regulation of tubulin expression (Fig 4A). Using RT-qPCR, we found that the levels of spliced TUBA1A and TUBB mRNAs (Fig 4B) were significantly higher in H1047R cells, suggesting post-transcriptional regulation of tubulin mRNA. In support of this result, transcription rate of TUBA1A and TUBB in H1047R cells was unchanged compared to the parental cell line (Fig 4C). To directly test whether PI3K was involved in the regulation of tubulin mRNA stability, we measured tubulin mRNA levels in parental and H1047R-mutant cells treated with the PI3K inhibitor GDC-0941 [28]. Upon PI3K inhibition in parental cell line, the stability of tubulin mRNAs remained unchanged (Fig 4B) despite the increased transcription rate (Fig 4C). Importantly, PI3K inhibition restored normal levels of spliced TUBA1A and TUBB mRNA in H1047R-mutant cell line, comparable to the parental cell line (Fig 4B), while it increased the transcription of TUBA1A, but not TUBB mRNA (Fig 4C). To test whether these results were more general, we measured TUBA1A and TUBB gene expression in another cell line carrying an activating mutation in the helical domain of PI3K (E545K [29]). We found that the levels of spliced TUBA1A and TUBB mRNA were significantly higher in E545K compared to parental cells (Fig 4B), while the rate of TUBA1A and TUBB transcription remained unchanged (Fig 4C), again consistent with regulation by mRNA stability. Treatment of E545K cells with GDC-0941 PI3K inhibitor restored normal levels of spliced TUBA1A and TUBB mRNA (Fig 4B), suggesting that PI3K activity regulated tubulin mRNA stability. As in parental and H1047R cells, transcription of TUBA1A and TUBB was higher in E545K mutant cells treated with GDC-0941 compared to control-treated cells, suggesting transcriptional activation triggered by PI3K inhibition (Fig 4C). We conclude that, despite activating transcription, PI3K inhibition post-transcriptionally suppresses tubulin mRNA levels, presumably through autoregulation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. PI3K activity increases tubulin gene expression through autoregulation.

(A) DGE in parental and cells expressing constitutively active PI3K H1047R mutant cells, treated with DMSO or indicated PI3K inhibitors for 4 h. (B and C) Relative expression of TUBA1A and TUBB unspliced pre-mRNA (B) and spliced TUBA1A and TUBB mRNA (C) in parental and PI3K mutant cell lines H1047R and E545K, treated with DMSO control or 1 μM GDC-0941 for 4 h. All the relative gene-expression data are normalized to housekeeping gene GAPDH or RPL19 and to DMSO-treated parental cells. (D) Tubulin partitioning to unpolymerized (soluble, S) and polymerized (P) normalized to loading control (GAPDH for S; HISH3 for P) and to DMSO-treated parental cells in parental and cells expressing constitutively active PI3K mutants H1047R and E545K, treated with DMSO or 1 μM GDC-0945 for 4 h. (E) DGE in A2058 cells treated with DMSO control or the indicated concentrations of PI3K inhibitors for 6 h (GEO series GSE66343 [28]). (F–G) RT-qPCR in A2058 cells treated with 1 μM GDC-0941 for indicated periods of time (x-axes). (H) Tubulin partitioning to unpolymerized (S) and polymerized (P) normalized to loading control (GAPDH for S; HISH3 for P) and to DMSO-treated cells in control and cells treated for 6 h with indicated 1 μM PI3K inhibitors. (I) Representative immunofluorescence images of control and cells treated with 1 μM GDC-0941 for 6 h and stained with anti-EB1 antibody (cyan), anti-β-actin, and Hoechst (red). (J) Relative number of detected EB1 comets per cell area in DMSO control and A2058 cells treated with 1 μM GDC-0941 for 6 h, normalized to DMSO-treated cells (100–200 cells). (K) A proposed model describes the mechanism through which PI3K signaling regulates tubulin gene expression through modification of tubulin autoregulation. Bar plots in all panels represent average values, and error bars standard deviations from three independent biological replicates. *p < 0.05, **p < 0.01, ***p < 0.001 in paired Student t test, in all panels. Black and white asterisk represent statistical significance relative to DMSO-treated parental cell line, green asterisks represent statistical significance relative to DMSO-treated control of the same cell line. a.u., arbitrary unit; DGE, differential gene expression; EB1, end-binding protein 1; GAPDH, Glyceraldehyde 3-phosphate dehydrogenase; GEO, Gene Expression Omnibus; HISH3, histone H3; pAKT, phospho-AKT; PI3K, phosphatidylinositol-4,5-biphosphate 3-kinase; RPL19, ribosomal protein L19; RT-qPCR, reverse-transcription quantitative PCR; TUBA, α-tubulin; TUBB, β-tubulin; TUBD, δ-tubulin; TUBE, ε-tubulin; TUBG, γ-tubulin; WT, wild type.

https://doi.org/10.1371/journal.pbio.3000225.g004

We next sought to determine whether PI3K activity regulated the stability of tubulin mRNAs via changes in microtubule stability and altered ratio of soluble tubulin dimer and polymer. To test whether PI3K activation causes a decrease in the concentration of soluble tubulin dimer, we performed biochemical tubulin partitioning on parental, H1047R, and E545K cells. We found that both H1047R and E545K cell lines had less soluble tubulin dimer (Fig 4D; note the differences in fractions “S” across the different cell lines) and more microtubule polymer than the parental cell line (Fig 4D; note the differences in fractions “P” across the different cell lines). PI3K inhibition in these cells increased levels of soluble tubulin dimer and decreased microtubule polymer, comparable to the levels observed in parental cells (Fig 4D; note the differences in fractions “S” and “P” between cells treated with DMSO and GDC-0941). These data are consistent with a model in which PI3K activity increases microtubule stability, leading to a decrease in soluble tubulin and an increase in tubulin mRNA stability via the autoregulation pathway. PI3K inhibition with a small molecule reverses these changes. While these data do not reveal the mechanism by which PI3K regulates microtubule stability, they are consistent with published models [30].

To generalize our findings, we investigated a second cell line, A2058, in which our bioinformatic analysis revealed down-regulation of tubulin mRNAs upon PI3K inhibition. PI3K is constitutively active in these cells (Fig 4H). We treated A2058 cells with the PI3K inhibitors BKM-120 [31] (inhibitor for p110α/β/δ/γ), BEZ-235 [32] (a dual PI3K and mammalian Target Of Rapamycin (mTOR) inhibitor for p110α/γ/δ/β), or GDC-0941 [33] (inhibitor for PI3Kα/δ). In each case, drug treatment decreased the levels of tubulin mRNA (Fig 4E). Performing RT-qPCR on cells treated with 1 μM GDC-0941 (IC90 [28]), we found that PI3K inhibition post-transcriptionally down-regulated spliced tubulin mRNAs (Fig 4F) while leaving intact the transcription rate of tubulin genes (Fig 4G), consistent with regulation through mRNA stability. Biochemical partitioning of tubulin revealed that treatment with 1 μM BKM-120, BEZ-235, or GDC-0941 indeed caused an increase in levels of soluble tubulin dimer and a reduction of microtubule polymer (Fig 4H). To further confirm microtubule destabilization upon PI3K inhibition, we used immunofluorescence-based microscopy. We fixed and stained control and PI3K-inhibited cells using anti-EB1 antibody and counted the number of EB1-positive microtubule plus-tips per cell area as a readout for growing microtubules. Our data showed significant reduction in the number of growing microtubules in cells treated with 1 μM GDC-0941 compared to control-treated cells (Fig 4I and 4J).

Several studies have identified off-target destabilization of microtubules by kinase inhibitors, so it was important to test for this possible artifact. In RPE1 hTert cells, with low basal PI3K activity, we did not observe changes in the number of EB1-positive microtubule plus-tips upon PI3K inhibition with GDC-0941 or BEZ-235, suggesting the absence of off-target activity on tubulin (S4A and S4B Fig). Consistent with previous reports of off-target effects on microtubules [34], high dose of PI3K inhibitor BKM-120 caused a reduction in the number of EB1-positive microtubule plus-tips (S4A and S4B Fig).

Taken together, we conclude that PI3K activity positively regulates tubulin levels in cancer cells via microtubule stabilization, which lowers the concentration of soluble tubulin dimer, thus increasing tubulin mRNA stability through the autoregulation pathway. Our findings show for the first time, to our knowledge, that tubulin autoregulation via mRNA stability mediates changes in total tubulin concentration triggered by an important signaling pathway. PI3K inhibition in cell lines with low constitutive kinase activity did not perturb microtubules (Fig 4D, S4A and S4B Fig), suggesting that this regulation is selective for cancer cells and normal physiological states in which PI3K is strongly activated.

Autoregulation suggests functions of tubulin beyond building microtubules

A striking feature of our DGE data in quiescent RPE1 hTert cells, which was also seen in two large GEO studies that we reanalyzed, is a lack of coregulation of other microtubule components with tubulin: tubulin autoregulation signature extends to all expressed tubulin but not to microtubule-associated, plus-tip–binding, or motor proteins, nor are these components coregulated on the transcriptional level in response to microtubule damage. This finding is quite different from coordinated gene expression in other biological responses; for example, many of the genes required to build lysosomes [38] or mitochondria [39] are coordinately regulated. The specific regulation of all expressed tubulins, but not other microtubule proteins, suggests that cells care more about the concentration of tubulin than other components of the microtubule cytoskeleton. Perhaps this simply reflects the central role of tubulin in building microtubules. A related puzzle is that autoregulation is counterhomeostatic from the perspective of microtubule mass, at least in response to microtubule drug perturbation. CA4 treatment reduces microtubule mass, yet cells respond by destabilizing tubulin mRNA and reducing protein levels, which would seem to exacerbate the problem, not correct it. The converse is true for PTX. An intriguing possible explanation for both the apparent paradoxes is that soluble tubulin dimer has some function other than building microtubules that necessitates tight control of its concentration. For example, soluble tubulin was proposed to gate Voltage-Dependent Anion Channels (VDACs) in mitochondria [40], and it is conceivable that autoregulation evolved to regulate this function, not microtubule assembly. Another possibility is that autoregulation evolved to balance the synthesis of TUBA and TUBB subunits, in which case other components of the microtubule cytoskeleton are irrelevant. But if this is the case, why is TUBG1 also coregulated? These are speculations, but they point out that we truly do not understand the adaptive benefit of tubulin autoregulation despite its frequent occurrence in cellular physiology, as revealed by our bioinformatics analysis.

Discovery of physiological processes that change microtubule stability from transcriptomic data

Research in the microtubule field has been driven largely by microscopy and biochemistry, while transcriptomic approaches have been neglected. Our CLIC search shows that bioinformatic mining of public data sets can be used as a starting point for discovery. As a positive control, this approach returned multiple studies of microtubule-perturbing agents, such as PTX and bisphenol-A. Unexpectedly, we uncovered many more previously unreported conditions that perturb tubulin gene expression that remain to be investigated.

Our mechanistic follow-up study showed that PI3K activity promotes stabilization of tubulin mRNA. This effect is achieved through the microtubule-stabilizing activity of PI3K signaling and tubulin autoregulation. Notably, the effect of autoregulation on basal tubulin levels, and on response to PI3K inhibitors, was stronger in cells carrying an activating mutation in PI3K. The inhibition of PI3K in wild-type and H1047R and E545K mutant MCF10A cells dramatically increased the transcription of TUBA1A and TUBB through mechanisms that remain to be explored. Suppression of TUBA1A- and TUBB-spliced mRNA stability through tubulin autoregulation counteracts the increased transcription rate, restoring normal levels of tubulin biosynthesis. Activating mutations in PI3K cause many changes in growth and metabolism in cancer cells [41]. Our data add change in tubulin levels to that picture, which could be significant for successful mitosis, response to microtubule drugs, and metastasis pathways regulated by microtubules.

Our genome-scale data analysis and PI3K investigation point to the importance of tubulin autoregulation in normal biology and cancer pathophysiology. Tubulin autoregulation was a central concern in the microtubule field in the 1980s, but it was abandoned as the field moved towards biophysical directions in the 2000s. We feel it is now important to investigate the full molecular basis of the autoregulation pathway and extend the pioneering work on the role of the N-terminal peptide of nascent TUBB [15].

Cell culture and drug treatments

All the cell lines used in this study were grown at 37°C with 5% CO₂ in a humidified incubator. hTert-RPE1-eGFP-EB3 (a kind gift from D. Pellman, HMS, Boston, MA, USA), and hTert-RPE1 (ATCC, Manassas, VA, USA) cell lines were grown in Dulbecco's modified medium (nutrient mixture F12, DMEM/F12) supplemented with 10% FBS and 1% (vol/vol) penicillin/streptomycin (pen/strep). Parental and PI3 kinase mutant MCF10A cell lines, H1047R and E545K (a gift from J. Brugge, HMS), were grown as previously described [29]. Prior to drug treatment, parental and PI3K mutant MCF10A cell lines, H1047R and E545K, were washed 2 times with 1× PBS and grown in medium lacking EGR [29]. Cells were incubated overnight. A2058 cells (a kind gift from P. Sorger, HMS), U2OS cells, and A549 cells (ATCC) were grown in DMEM supplemented with 10% FBS and 1% (vol/vol) pen/strep. Microtubule drugs were used at 1 nM or 100 nM for CA4 or 3 nM or 300 nM for PTX for 6 and 24 h. PI3K inhibitors BKM-120, BEZ-235, and GDC-0941 (a gift from Nathanael Gray, HMS) were administered at 1 μM. All the microtubule poisons and PI3K inhibitors were dissolved in DMSO, and 0.01% DMSO was used as vehicle control. D-glucose–free medium was prepared using D-glucose–free DMEM (Thermo Fisher Scientific, Waltham, MA, USA) supplemented with 10% dialyzed FBS (Thermo Fisher Scientific), and 1% vol/vol pen/strep (Thermo Fisher Scientific). Media with D-glucose were prepared by adding D-glucose (Sigma Aldrich, St. Louis, MO, USA) to the glucose-free medium at 4.5 g/l final concentration. L-glutamine–free medium was prepared using L-glutamine–and sodium-pyruvate–free DMEM (Corning, Corning, NY, USA) supplemented with 10% dialyzed FBS, 1% sodium pyruvate (Corning), and 1% vol/vol pen/strep. Media with L-glutamine were prepared by adding L-glutamine (Corning) to the L-glutamine–free medium at 2 mM final concentration. For D-glucose and L-glutamine deprivation, standard growth medium was removed, cells were washed 3 times with 1× PBS, and media with (control) or without D-glucose or L-glutamine were added to cells.

Flow cytometry cell-cycle profiling

In the flow cytometry validation experiments, RPE1 hTert cells were grown per treatment condition on 10-cm petri dishes for 1 day for cycling culture and 5 days after reaching full confluency for quiescent cultures. Cells were then treated with DMSO control, 100 nM CA4, or 300 nM PTX for 6 h. For each condition, cells were detached from the dish in 0.05% Trypsin solution, resuspended in 1× PBS supplemented with 1% FBS and spun down by centrifugation (400 × g for 5 min), followed by dispersion of the pellet into a single-cell suspension in 1× PBS. Cells were then fixed in 4% PFA for 20 min on ice, stained with DAPI (#D1306; Invitrogen, Carlsbad, CA, USA), and analyzed on a BD LSRII flow cytometer (BD Biosciences, Franklin Lakes, NJ, USA). The dye was excited with a 405-nm laser, and emitted fluorescence was detected with a 450/50 bandpass filter. For data analysis, FlowJo v.X.0.7 and the built-in cell-cycle quantification platform were used (the univariate model without any adjustments).

RNA-seq and data analysis

Total RNA was collected and purified using the PureLink RNA Mini Kit (Invitrogen, Thermo Fisher Scientific). RNA concentration and quality were determined using NanoDrop and Bioanalyzer, respectively, and 500 ng of purified RNA was used as input for the Illumina TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA). Barcoded libraries were pooled and sequenced by The Bauer Core Facility at Harvard University, where the pooled library was quantitated using KAPA and single-end sequenced on an Illumina NextSeq (Illumina). RNA-seq reads were mapped using STAR [42] (version 2.1.0j) and processed using HTSeq-count [43] (version 0.6.1). GRCh38 reference genome and transcript annotations were used for gene mapping; Entrez Gene identifiers and the org.Hs.eg.db database were used for genome-wide annotation [44]. DGE and statistical analysis were performed using the edgeR package [45]. Genes with >50 CPM and a fold change significantly different from zero in Wilcoxon signed-rank test (p < 0.05) were marked as differentially expressed genes, based on three biological replicates. EnrichR was used for GSEA [46,47], and gene sets with Benjamini–Hochberg adjusted p-values p < 0.01 were considered statistically significant.

Microarray data analysis

Publicly available microarray data were read, normalized, and analyzed using linear models and empirical Bayes methods for assessing differential expression (limma, [48]). The microarray probes that hybridize with tubulin genes were extracted, and heatmaps of differential expression were plotted using R.

Bioinformatic analysis

We used the CLIC-gene online tool for bioinformatic analysis of tubulin expression correlation [26]. Querying a list of all tubulin genes against human and mouse data set platforms, the algorithm calculated average Pearson expression correlation for each data set within one platform separately and using expression correlation coefficients obtained from the top-ranked 417 human and 122 mouse DGE data sets. The 417 human data sets with high tubulin expression correlation were then manually annotated.

RT-qPCR

Reverse transcription was performed from 500 ng of purified total RNA, using SuperScript IV (Invitrogen) and random hexamer primers, according to the manufacturer’s protocol. RT-qPCR reaction was performed using 5 ng of cDNA and 2× SYBRGreen master mix (Thermo Fisher Scientific) on a BioRad thermocycler (BioRad, Hercules, CA, USA). For each reference and gene of interest, two sets of primers were designed: one set of primers amplified specifically unspliced pre-mRNA, while the other set of primers amplified specifically spliced mRNA. Primer sequences are listed in S1 Table. All primer pairs were validated by PCR, followed by gel electrophoresis. PCR products were further submitted to Sanger sequencing and subsequently mapped against the genome, confirming correct product amplification. RT-qPCR data analysis was performed using the ddCt method [49]. Statistical significance was determined using two-tailed paired Student t test and p < 0.05.

Western blotting and tubulin partitioning

Whole-cell extracts for immunoblot analysis were prepared by cell lysis in 2× SDS sample buffer containing 50 mM Tris-HCl (pH 6.8), 2% SDS, 10% glycerol, 1% β-mercaptoethanol, 12.5 mM EDTA, and 0.02% bromophenol blue in water. Samples were denatured at 100°C for 10 min prior to SDS-PAGE gel electrophoresis on 1.5 mm NuPAGE Novex 4%–12% Bis-Tris Protein Gels (Thermo Fisher Scientific), using the Mini Trans-Blot Cell system (BioRad).

For tubulin partitioning, cells were grown in 12-well dishes to 70%–80% confluence. Prior to lysis, cells were washed 1× with 1× PBS prewarmed to 37°C, and PBS was aspirated. Unpolymerized tubulin was extracted in 300 μl of tubulin extraction buffer containing 60 mM PIPES (pH 6.8), 25 mM HEPES (pH 7.2), 10 mM EGTA (pH 7–8), 2 mM MgCl₂, 0.5% TritonX-100, 10 μM PTX (Sigma Aldrich), and a protease inhibitor tablet (Roche, Branchburg, NJ, USA) for 2 min at room temperature. Extraction buffer with solubilized protein was collected and mixed with 100 μl 4× SDS sample buffer. The remaining material, containing polymerized tubulin, was then lysed with 400 μl 1× SDS sample buffer. Samples were denatured at 100°C for 10 min. Prior to SDS-PAGE gel electrophoresis, equal aliquots of samples were concentrated 2× by evaporation.

For all immunoblots, Precision Plus Protein Dual Color Standard ladder was used (BioRad). Dilutions and primary antibodies used: 1:15,000 GAPDH (#2118; Cell Signaling Technology, Danvers, MA, USA), 1:10,000 TUBA (DM1alpha, #05–829; MilliporeSigma, Burlington, MA, USA), 1:10,000 Histone H3 (#4499; Cell Signaling Technology), and 1:500 phospho-Akt S437 (#9271; Cell Signaling Technology). The following secondary antibodies were used: 1:15,000 conjugated goat anti-mouse DyLight 680 (#35518; Thermo Fisher Scientific) and 1:15,000 conjugated goat anti-rabbit DyLight 800 (#35571; Thermo Fisher Scientific). All the blots were imaged and analyzed using The Odyssey Infra-Red Imaging System (LI-COR, Lincoln, NE, USA) equipped with Image Studio software. Statistical analysis and data plotting were performed in Microsoft Excel.

Quantitative immunofluorescence

For quantitative immunofluorescence, cells were seeded and grown on coverslips overnight, prior to drug treatment. For immunostaining of EB1, cells were washed once with 1× PBS at 37°C and fixed in cold methanol at ^l20°C for 20 min. Cells were then permeabilized in 4% paraformaldehyde in PBS at room temperature for 20 min. After blocking in 2% bovine serum albumin in PBS, primary antibodies were incubated at room temperatures for 1 h. The following primary antibodies and dilutions were used: 1:500 mouse anti-EB1 (#2164; Cell SignalingTechnology), 1:1,000 rabbit anti-β-actin (#4970; Cell Signaling Technology), 1:15,000 Hoechst 3342 (#H3570; Thermo Fisher Scientific). Cells were then stained with secondary antibodies at room temperature for 30 min. The following secondary antibodies and dilutions were used: 1:400 Alexa Fluor goat anti-rabbit IgG (#11008; Life Technologies, Carlsbad, CA, USA), and 1:400 Alexa Fluor goat anti-mouse IgG (#11001; Life Technologies). Three-dimensional image stacks of interphase cells were acquired in 0.2-μm steps using a 60× NA 1.42 objective on an Olympus DeltaVision microscope (GE Healthcare, Chicago, IL, USA) equipped with a DAPI/FITC/TRITC/CY5 filter set (Chroma, Bellows Falls, VT, USA), and an sCMOS 5.5 camera (PCO, Kelheim, Germany). EB1 comet counting and cell area measurements were performed on maximum-intensity z-stack projections using a custom-made ImageJ macro, which is available upon request. Statistical analysis and data plotting were performed in R.

Statistical analysis

For gene-expression profiling, statistical analysis was supplied as part of the edgeR package and performed according to the user manual. For tubulin gene-expression correlation, statistical analysis was performed using the online CLIC-gene tool. All the Student t test analyses were performed considering two-tailed distribution and two degrees of freedom.