Gene length as a biological timer to establish temporal transcriptional regulation

ABSTRACT

Transcriptional timing is inherently influenced by gene length, thus providing a mechanism for temporal regulation of gene expression. While gene size has been shown to be important for the expression timing of specific genes during early development, whether it plays a role in the timing of other global gene expression programs has not been extensively explored. Here, we investigate the role of gene length during the early transcriptional response of human fibroblasts to serum stimulation. Using the nascent sequencing techniques Bru-seq and BruUV-seq, we identified immediate genome-wide transcriptional changes following serum stimulation that were linked to rapid activation of enhancer elements. We identified 873 significantly induced and 209 significantly repressed genes. Variations in gene size allowed for a large group of genes to be simultaneously activated but produce full-length RNAs at different times. The median length of the group of serum-induced genes was significantly larger than the median length of all expressed genes, housekeeping genes, and serum-repressed genes. These gene length relationships were also observed in corresponding mouse orthologs, suggesting that relative gene size is evolutionarily conserved. The sizes of transcription factor and microRNA genes immediately induced after serum stimulation varied dramatically, setting up a cascade mechanism for temporal expression arising from a single activation event. The retention and expansion of large intronic sequences during evolution have likely played important roles in fine-tuning the temporal expression of target genes in various cellular response programs.

KEYWORDS:

• enhancer elements
• evolution
• intron size
• nascent RNA sequencing
• serum response
• transcription initiation

Introduction

Human genes come in a wide variety of lengths; protein-coding genes range in length from less than a hundred to over 2 million base pairs.¹ One consequence of this extensive range in gene size is that the time needed to complete transcription elongation of a full-length transcript is also highly variable among genes. Hypothetically, at a constant transcription elongation rate of 1.4 kb/min,^2-5 a 100 bp gene would take only a few seconds to finish while a 2 Mb long gene at the same rate would take over 24 hours to complete elongation. These gene size differences are largely driven by variations in the lengths of intronic sequences. Long genes tend to be mostly comprised of intronic sequences that are transcribed but then excised by the spliceasome and degraded by the nuclear RNA exosome.⁶ Additionally, introns tend to be longer in humans than in other vertebrate species,⁷ and longer in tissue-specific genes compared to housekeeping genes.⁸ Furthermore, very long genes promote genomic instability through increased risks for the formation of copy number variations (CNVs) and common fragile sites.⁹ Why have long genes survived evolutionary pressures for shortening despite the time and energy burden and the threat to genomic integrity they present? While introns undoubtedly play critical roles in genome organization, hosting ncRNA genes, RNA processing and isoform diversity,^10-12 their contribution to gene size variability suggests that intron length may also serve an important biological function in regulating the temporal expression of genes induced simultaneously in response to stress or cell stimuli.

Early observations of large introns in Drosophila developmental genes led to the “intron delay hypothesis” which postulates that intron length may play a role in gene expression timing.¹³ Indeed, it has been demonstrated that during Drosophila early development, the gap gene kni is expressed while its cognate gene knrl, which has more intronic sequence, does not have enough time to be expressed due to rapid cell divisions disrupting its transcription.¹⁴ More recently, it was shown that intronic length is important during somite segmentation of mouse embryos.¹⁵ During this time, the Hes7 gene is cyclically expressed as a result of regulation through a negative feedback loop involving the Hes7 protein. Removal of introns from Hes7 led to earlier expression and abolishment of oscillatory expression, and resulted in developmental defects. While intron length is implicated in the regulation of expression timing for certain genes during development, whether variable intron lengths are important for coordinating gene expression programs outside of development is unknown.

Gene expression programs allow cells to react to specific changes in their environment through temporally coordinated transcription of response genes. Various signaling cascades activate preexisting transcription factors, which in turn induce primary response genes. Immediate-early genes induced as part of the primary response are often rapidly and transiently expressed, and many of these genes encode transcription factors.¹⁶ The target genes of these transcription factors are part of the secondary response in which gene induction requires de novo protein synthesis and therefore necessitates a delay before transcription initiation of these genes. Initiation of transcription is classically viewed as a key regulatory point for the timing of gene expression, however other transcriptional events such as pre-positioned RNA polymerase II (RNAPII) at promoter proximal pause sites and rate of transcriptional elongation which is affected by intron base composition, repetitive sequences, and RNA splicing complex formation, can also influence temporal expression patterns during the early activation stages of gene programs.^17,18

The early response to serum stimulation provides an excellent model for studying immediate-early gene expression programs. While this response involves many triggers and signal transduction pathways making it difficult to elucidate mechanisms responsible for transcriptional regulation, it is a well-studied response for inducing robust transcriptional reprogramming. Fibroblasts grown in serum-free media enter a G₀ state of minimal metabolic activity, and subsequent addition of serum to these cells results in global gene expression changes as cells prepare to re-enter the cell cycle.¹⁹ While the genome-wide early serum response has been previously explored, these studies focused on expression changes based on steady-state mRNA levels.^20,21 Because such analyses only provide information about the steady-state of completed transcripts, they do not directly indicate timing of upstream transcriptional events such as induction or repression of transcription initiation. Some studies have explored nascent transcription of early response genes using qRT-PCR,^17,22 but this approach is limited to the assessment of transcription of a handful of selected genes. We recently used nascent RNA Bru-seq to assess the dynamics of transcription initiation and elongation during the first 2 hours after serum stimulation.²³ In the present study we used Bru-seq and BruUV-seq to assess induction and repression of transcription initiation, and changes in the activity of enhancer elements immediately after serum stimulation.^24-26 Because Bru-seq allows for the assessment of transcription initiation without the need for the generation of a mature mRNA product, our study provides a comprehensive profile of the immediate transcriptional response to serum. Our results confirm the rapid induction of many previously described serum-response genes and provide a novel list of genes immediately repressed by serum stimulation. We also show that many of the highly induced genes were associated with activated putative enhancers while downregulated genes were rarely associated with inactivated enhancers, suggesting that enhancer regulation primarily plays a key role in the regulation of rapid gene induction. Finally, serum-induced genes were found to be on average significantly longer than housekeeping genes and serum-repressed genes, highlighting the importance of gene size variability and the role of gene length in temporal global expression timing.

Results

Identification of immediate serum response genes using Bru-seq

To investigate the immediate transcriptional response to serum, we compared nascent RNA expression in starved and serum activated normal human fibroblasts using Bru-seq.^23-25 Cells were grown in serum-free media for 48 hours. Serum was then added back to the media (or not for the starved control), and nascent RNA was immediately labeled with bromouridine (Bru) for 30 minutes (Fig. 1A). The Bru-RNA was isolated, used to prepare cDNA libraries, deep sequenced, and mapped to the reference genome (hg19). The experiment was repeated and the transcription profiles of the 2 biological experiments were highly correlated (Fig. S1). The mapped reads typically give a jagged appearance, which is quite reproducible and most likely stem from an inherent feature of library preparations primarily caused by variations in PCR amplification. Among the genes that were immediately induced, we observed known immediate-early genes such as the transcription factor gene FOS (Fig. 1B). This was evident through an increase in reads across the entire gene. In long induced genes, such as PDE7B (∼344 kb) which encodes a nucleotide phosphodiesterase, the increase in reads was present at the beginning of the gene but only extended partially into the body of the gene (Fig. 1C). We previously measured a median elongation rate of ∼1.4kb/min in the first 40 kb of genes of our fibroblast cell line.² Consistent with this rate, the 30-minute labeling period did not allow enough time for newly initiated PDE7B transcripts to be completed. In contrast to RNA techniques based on steady-state RNA analysis, the Bru-seq approach allowed us to identify genes that were transcriptionally repressed immediately after serum addition, such as TRIB2 (Tribbles pseudokinase 2) (Fig. 1D). For long repressed genes, such as the signaling gene GNG2 (guanine nucleotide binding protein 2) (109 kb), we observed a decrease in reads at the beginning of the gene and a receding wave of reads toward the 3′ end of the gene (Fig. 1E). This receding wave corresponds to transcripts initiated prior to serum addition that continued to elongate toward the 3′-end of the gene during the labeling period.

Figure 1. (A) Experimental outline. Bromouridine (2 mM) labeling was performed for 30 minutes on serum starved human fibroblasts or together with serum for 30 min on previously starved human fibroblasts. Nascent RNA sequencing reads expressed as RPKM are shown for (B) FOS, (C) PDE7B, (D) TRIB2 and (E) GNG2 in starved cells (orange trace) and in cells serum stimulated for 30 min (blue trace).

For our genome-wide analysis, we calculated nascent RNA expression levels for all genes in starved cells and during the 30-minute serum stimulation period. For genes over 30 kb, we calculated expression levels based on the first 30 kb of the genes because we predicted that newly induced or repressed long genes would exhibit read changes within this region but not at the end of these genes. Out of the 6958 genes fitting our expression criteria (at least 300 bp long and expressed above 0.5 RPKM), we identified 873 significantly induced and 209 significantly repressed genes (adjusted p-values < 0.05, n = 2) during the first 30 min of serum stimulation (Table S1). Our results demonstrate that serum stimulation of starved human fibroblasts results in extensive, immediate transcriptional changes that are detectable by Bru-seq, including a novel set of immediately repressed genes.

Comparisons between Bru-seq and microarray technologies for serum-induced transcription timing

Previously defined groups of immediate-early genes based on increases in steady-state mRNA levels tend to be dominated by small genes.^16,22 This is partially due to timing constraints on elongation since these techniques are biased toward detection of full-length transcripts. Bru-seq allows for the assessment of instantaneous changes in transcription, and by analyzing reads just downstream of transcription start sites (TSSs), changes in transcription initiation can be inferred independently of gene length. Such immediate assessment of transcriptional changes is not afforded by microarray or conventional RNA-seq analysis because changes are slow to appear in steady state RNA pools due to the “noise” of pre-existing RNA. Furthermore, if RNA is captured through poly(A) selection, changes will only be detected after completion of full-length transcripts.

A previous study used the microarray approach and polyA+ enriched RNA to assess the timing of transcriptional changes in human fibroblasts following serum stimulation,²⁰ and we compared this dataset to our Bru-seq dataset. For short genes, both Bru-seq and the microarray assay detected rapid increases in transcription levels by the first 30 minutes following serum stimulation (Fig. S2). In contrast, for longer genes, the microarray analysis detected increased transcript levels at later labeling periods following serum addition while Bru-seq detected an immediate increase of nascent transcription at the beginning of the gene within the first 30 minutes. This timing discrepancy can be explained by the delay between productive initiation of RNAPII and the time needed to elongate across the entire gene to produce full length transcripts detected by the microarray assay in the total RNA pool. These comparisons of transcriptional response timing demonstrate how different methods capture distinct stages of transcription.

Gene size-dependent induction of functional pathways

We next tested the hypothesis that gene length may play a functional role in the temporal regulation of gene expression during cellular responses. Although a large set of genes may be induced simultaneously, variations in gene size allow for staggered generation of full-length RNA products. This would lead to temporally ordered translation of mRNAs at ribosomes and staggered production of proteins involved in different response pathways. We performed gene set enrichment analysis (GSEA) on our list of expressed genes which were pre-ranked according to fold change in the serum-stimulated sample compared to the starved sample (Table S2). The pathways which showed enrichment among genes induced by serum stimulation included “TNF signaling via the transcription factor NFκB,” “TGFβ signaling,” “apical junction,” and “hypoxia” (Fig. 2A). Pathways which showed enrichment among genes suppressed by serum included “meiotic recombination,” “RNA pol I transcription,” “peroxisome,” and “oxidative phosphorylation” (Fig. 2B). To examine whether induced genes of different sizes were enriched for specific biological pathways, we grouped the 873 serum-induced genes based on size and performed DAVID functional annotation analysis (Table S3, Figs. 2C-G, S3). The induced gene groups were compared to a background of all expressed genes within the same size category in order to eliminate biases of pathways that may be overrepresented in certain size categories. Pathways such as “focal adhesion” and “actin cytoskeleton” were enriched in all size categories while “MAPK signaling pathway” was enriched in 4 out of 5 size groups. “Regulation of transcription” was enriched in the 3 smallest size classes, while “pathways in cancer” was enriched in the 2 largest gene categories. For genes repressed by serum stimulation, we grouped the 209 genes into 2 size groups, smaller or larger than 15 kb. Among the smallest gene group, enriched gene sets included “transcription,” “mitochondrial oxidative phosphorylation,” and “chromatin assembly” (Fig. 2H), while for the larger genes, enrichment for “actin cytoskeleton,” “focal adhesion,” “pathways in cancer,” and “endocytosis” was found (Fig. 2I).

Figure 2. Gene set enrichment analysis of immediate serum-response genes. Normalized enrichment scores calculated by GSEA for pathways identified in induced (A) and repressed (B) serum response genes. Enriched pathways identified by DAVID in induced (C-G) and repressed (H-I) gene sets grouped based on gene size. The data is expressed as -log10 of p-values from two independent biological experiments.

Primary transcription factors responding to serum stimulation

To explore what transcription factors may be responsible for regulating the immediate serum response in human fibroblasts, we performed GSEA analysis to look for transcription factor-binding motifs enriched among the induced genes (Table S4). As expected, the transcription factor binding motif for serum response factor (SRF) had the highest enrichment score, followed by CREBP1 (ATF2), and other members of the ATF/CREB transcription factor family (Fig. 3A-C). Among repressed genes, we found enrichment of 3 motifs that have not yet been associated with any known transcription factor (Fig. 3D-F). Because our Bru-seq set of immediately repressed genes is novel, much less is known about the transcription factors and regulatory mechanisms controlling the suppression of these genes.

Figure 3. Transcription factor binding motif enrichment among immediate serum response genes. Normalized enrichment scores (NES) calculated by GSEA for highest scoring transcription factor binding sites identified in induced (A) and repressed (B) serum response genes.

Rapid activation of enhancer elements following serum stimulation

We recently developed BruUV-seq, which allows for the genome-wide identification of putative active enhancer elements.²⁶ Here we used BruUV-seq to investigate the relationship between immediately induced serum response genes and concurrently activated enhancer elements. In BruUV-seq, cells are exposed to UV-irradiation prior to Bru-labeling, resulting in increased signal for certain classes of non-coding RNAs such as RNA generated at active enhancer elements (eRNA). This enhanced eRNA signal is likely due to UV-induced transcription-blocking lesions and inhibition of RNA exosome components, which together allow us to capture unstable RNA species.^26-28 We identified several induced genes that had corresponding increases in activity at nearby putative enhancer elements. For example, Bru-seq identified a 9-fold induction of the FOS gene (Table S1), which was accompanied by an increase in reads at nearby UV-enriched peaks that overlap with an annotated lncRNA (Fig. 4A). These peaks were not seen in the Bru-seq data, reflecting the unstable nature of these lncRNAs, presumably due to rapid degradation by the RNA exosome, which suggests that they may represent eRNAs.²⁶ Similar findings of nearby putative enhancer elements were made for the highly induced genes NR4A1, TNC, ID1, ID2, and ID3 (Fig 4B-F). Out of the top 50 most highly induced genes, we found evidence that 39 of these were adjacent to at least one putative activated enhancer element, and out of these, 18 genes were adjacent to 4 or more putative activated enhancer peaks (Table S1). In sharp contrast, only one of the top 15 most highly repressed genes following serum stimulation coincided with suppression of nearby putative enhancer elements (Table S1). Thus, during the serum response, immediate induction of many genes appears to be associated with rapid activation of nearby enhancer elements, while genes rapidly repressed by serum stimulation appear to be inhibited in an enhancer-independent manner.

Figure 4. Putative enhancer activation following serum stimulation. Bru-seq (top) and BruUV-seq (bottom) traces are shown for starved (blue) and serum stimulated cells (orange) for (A) FOS, (B) NR4A1, (C) TNC, (D) ID1, (E) ID2, and (F) ID3. Genes are shown on top in green for plus strand genes and red for minus strand genes. Transcripts and histone mark peaks from ENCODE data for normal human lung fibroblasts (NHLF) are shown below. Red arrows point at intergenic transcription peaks enhanced in BruUV-seq compared to Bru-seq, which align with peaks for the enhancer marks H3K4me1 and H3K27ac.

Diverse sizes of serum-induced transcription factor and miRNA genes

A prominent group of immediately induced genes following serum stimulation identified by Bru-seq were transcription factor genes (Table S4). Of the 873 immediately induced genes, 111 were transcription factor genes, while 21 of the 209 repressed genes encoded transcription factors. These transcription factor genes were of various sizes and expected to complete transcription at different times. Based on size, small transcription factor genes such as FOSB, NR4A1, and NR4A2, are expected to produce full-length transcripts within about 6 minutes (Fig. 5A-C), medium-sized genes such as BACH1, TSC22D2, and ZNF131 are expected to produce full-length transcript after about 30 minutes (Fig. 5D-F), and large genes such as NFKB1, KLF7, and HIVEP2 are expected to produce finished transcripts after 1-3 hours following serum addition (Fig. 5G-I).

Figure 5. Transcriptional delays as a result of gene length. Bru-seq traces for (A) FOSB, (B) NR4A1, (C) NR4A2, (D) CDK7, (E) REL, (F) LIMA1, (G) PDLIM5, (H) BTAF1, and (I) UBR4 during starved conditions (orange trace) and following serum stimulation (blue trace), along with estimated times required for transcription completion (right) based on length and an elongation rate of 1.4 kb/min.

MicroRNAs (miRNAs) are small non-coding RNAs that can target specific mRNAs for degradation or inhibit translation.²⁹ Annotation of miRNA genes has been difficult because rapid processing of miRNA precursors into mature miRNA do not allow traditional RNA-seq techniques to capture miRNA primary transcripts. We have used Bru-seq to identify and annotate primary miRNA gene sequences genome-wide in multiple cell lines and found that they vary dramatically in size (manuscript in preparation). In this study, we found several miRNAs genes that were transcriptionally upregulated in response to serum in human fibroblasts. Examples of miRNA genes that were rapidly induced following serum stimulation are the gene pairs MIR143:1/ MIR143:2 (Fig. 6A) and MIR30B/MIR30D (Fig. 6B). Other induced miRNA genes appeared to be regulated by transcription run-on past the termination sites of the upstream genes, such as MIR21 located downstream of VMP1 (Fig. 6C) and MIR3193 located downstream of ID1 (Fig. 6D). Because these primary miRNA transcripts are short, they are expected to be completely transcribed within the 30 min labeling period. In contrast, a miRNA cluster on chromosome 14 consisting of 43 miRNAs is transcribed from a single promoter as one 200 kb unit,²⁶ and therefore these miRNAs are subject to a delay in production. These variably sized primary miRNA transcription units set up a type of cascade mechanism for temporal gene regulation similar to transcription factor genes. The size of the transcription unit will determine the length of time delay before these transcripts can be processed into mature miRNAs for inhibition of translation and degradation of target mRNAs.

Figure 6. Regulation of transcription of miRNA genes of different sizes by serum stimulation. (A-E) Bru-seq traces for various miRNAs during starved conditions (orange trace) and following serum addition (blue trace).

Relative gene length is evolutionarily conserved

One possible role for long introns is to delay the completion of full-length transcripts. The simultaneous activation or repression of a set of genes will result in different completion times determined by gene length. We reasoned that if the serum response utilizes a gene length delay mechanism to regulate temporal dynamics, serum-response genes will be, on average, longer than housekeeping genes that are constitutively transcribed. To test this, we compared the median lengths among different human gene groups. The median length of transcription factor genes expressed in human fibroblasts is ∼27 kb (488 genes), similar to the median length of all genes expressed in these cells, which is ∼29.5 kb (6958 genes) (Fig. 7A). In contrast, human housekeeping genes expressed in human fibroblasts (as defined by a study which identified as genes expressed across 19 normal tissues¹⁰) were found to be significantly shorter with a median length of just ∼9.4 kb (336 genes). Next, we compared the median length of the serum induced genes to serum repressed genes and found that induced genes were overall significantly longer than repressed genes, suggesting that temporal expression timing of induced genes relies on gene size more than repressed genes. This also held true when examining the subgroup of transcription factor genes responding to serum stimulation. We estimated transcriptional completion timing of serum response transcription factors based on length (Fig. 7B). While we estimate that the majority of the transcription factors complete transcription within the first hour after serum stimulation, there are several transcription factors (27/131) that are likely to complete transcription later, some not until 2 hours following serum addition. Following translation, induced transcription factors likely activate another set of target response genes, and expression timing of these target genes would depend on gene size as well. This sets up a complex regulatory mechanism by which successive series of gene induction events set up precise temporal expression patterns as a function of gene size. We found that expressed genes in human fibroblasts and corresponding mouse orthologs exhibited similar size distributions when grouped as above (Fig. 7C). Overall, relative gene size is well-conserved between human and mouse (Fig. 7D). Previous studies have shown that genes tend to be larger in humans compared to other animals, largely due to intron lengthening.^7,30,31 A study comparing orthologous introns between human and mouse found that 70% of human introns are longer than corresponding mouse introns.³⁰ However, although human introns are longer, there is a strong conservation of relative lengths of specific introns. Therefore, the longest intron of a gene in human tends to be the longest intron of a gene in mouse. These results show that gene size distribution is conserved and suggest that the maintenance of relative gene size may be important for proper expression timing.

Figure 7. Serum response genes and their mouse orthologs exhibit a broad range of gene sizes. (A) Boxplots displaying the distributions of gene sizes for various gene sets: all genes expressed in either starved or serum stimulated cells (6958 genes), transcription factor genes expressed in either starved or serum stimulated cells (488 genes), housekeeping genes expressed in either starved or serum stimulated cells (336 genes), genes induced after serum stimulation (873 genes), genes repressed after serum stimulation (210 genes), transcription factor genes induced after serum stimulation (111 genes), and transcription factor genes repressed after serum stimulation (21 genes). Stars indicate significance of Mann-Whitney-Wilcoxon statistical tests comparing the medians of the different groups, with p-values *< 2 × 10^-5, **< 2 × 10^-7 and ***< 3 × 10^-16. (B) Induced and repressed transcription factor genes displayed according to estimated time needed to complete transcription. (C) Boxplots displaying the distributions of gene sizes for mouse orthologs of the gene sets shown in A. Only orthologs with one-to-one homology were assessed. All genes (5719 genes), all TFs (385 genes), housekeeping genes expressed (272 genes), induced genes (753 genes), repressed genes (169 genes), induced TFs (102 genes), and repressed TFs (17 genes). (D) Correlation plot of sizes of all human protein-coding genes and sizes of their mouse orthologs (15845 genes). Only orthologs with one-to-one homology were assessed.

Temporal expression for optimal protein complex assembly

Assembly of multi-protein complexes may benefit from temporal synthesis of individual components. To explore this, we focused on the GO term category “cellular components,” which includes complexes that we found to be enriched among our induced genes such as “actin cytoskeleton” and “plasma membrane” (Fig. 2C-F). Many of these complexes contained genes within several size categories (Fig. S3). Examples of transcriptional induction of 4 genes of varying lengths that encode actin cytoskeleton proteins are shown in Fig. 8. The ACTC1 gene contains short introns allowing for rapid transcription after induction, while the other 3 genes contain large introns causing transcriptional delays. These genes range in size from 7.6 kb for ACTC1 (Fig. 8A) to 122 kb for VCL (Fig. 8D), which translates into an estimated transcript completion times of 5 and 81 minutes, respectively. String interaction networks for the 4 different gene products are provided to show their physiological interactions.³² It would be of interest to study the relationship between gene expression timing and the temporal assembly of multi-protein complexes, such as actin cytoskeleton. Perhaps the temporal gene expression pattern driven by evolutionary selection of particular gene sizes could teach us about the importance of the order by which these complexes are assembled. Additional studies are needed to elucidate the importance of temporal expression of different protein components for optimal assembly of actin cytoskeleton and other structural complexes in cells entering the cell cycle following serum stimulation or exposed to other environmental cues.

Figure 8. Transcription of actin cytoskeleton genes following serum stimulation. Bru-seq traces for (A) ACTC1, (B) DSTN, (C) ACTN4, and (D) VCL during starved conditions (orange) and following serum addition (blue). String interaction networks are shown below.

Figure 9. Model for the role of gene length in establishing temporal expression patterns following serum stimulation. (A) Transcription induction occurs simultaneously for several genes (blue arrow), but gene length influences the completion timing of the transcript (red lines). These immediately induced transcription factors (colored circles) then go on to activate their own gene targets (red arrows), whose expression timing is also influenced by gene length. (B) Temporal, staggered expression timing established by gene length.

Discussion

Initiation of transcription is a key regulatory step for gene expression, where recruitment of RNAPII and promoter proximal pausing contribute to variations in temporal expression patterns during early cellular responses.^17,18 Gene length and transcription elongation rates also act as regulatory mechanisms controlling gene expression timing. Human genes span in size from a few hundred to over 2 million base pairs. The presence of introns in human genes necessitates splicing, which allows for the generation of alternative splice variants leading to the generation of multiple protein forms encoded by an individual gene. This diversification of the proteome is an important function of introns, but why do introns exhibit such diverse sizes? The time and energy exerted while transcribing large introns as well as the increased threat of genomic instability from very large transcription units⁹ would be good reasons for evolutionary pressures to shape genes into more compact sizes. Instead, the retention and expansion of large and variably-sized introns within genes could have been selected to stagger the production of their full-length transcripts following transcription initiation to temporally regulate specific gene programs.

Immediate-early response genes have been typified by their short length and rapid induction. Our genome-wide profiling of nascent transcription using Bru-seq indicated that for the serum response, these common traits are not due to selective initiation of small genes but instead are a result of elongation constraints on longer genes which results in the staggered generation of full-sized transcripts. If genes in a particular response pathway were of similar sizes, the completion of their mRNAs would be expected to occur at the same time, and this could perhaps cause a “traffic jam” at ribosomes. Instead, evolutionary selection of different sized genes within a specific response pathway allows for temporal spacing of mature mRNA arrival at ribosomes for protein synthesis. Considering that many immediate response genes encode transcription factors that in turn regulate sets of different-sized genes, this potentially allows for one initial triggering event to induce a complex network of temporal expression patterns through combinatory regulation by gene length and signaling cascades (Fig. 9). Additionally, induction of primary miRNAs transcription units of different sizes is expected to temporally space the production of mature miRNA and the subsequent downregulation of target gene expression. We predict that immediately transcribed genes for other cellular responses would also have a broad size range based on the complexity of the induced pathways and whether the response occurs over a longer period of time. In our previous study, we treated cells with TNF for 1 hour to stimulate an immune response,²⁴ and induced genes identified by Bru-seq in that study were also longer overall compared to all expressed genes (data not shown).

Previous studies investigating gene expression changes following serum stimulation have been performed using steady-state RNA, and in many studies the mature RNA was obtained via poly(A) tail selection. One advantage of using nascent RNA Bru-seq for gene expression analysis is that the data is not confounded by the presence of pre-existing RNA. As a result, the Bru-seq analysis captured both immediate induction and repression of transcription. The data obtained for immediately repressed genes is novel because it reflects transcription rate changes and is not linked to the turnover of pre-existing RNA. We identified certain gene ontology terms not previously linked to repression following serum addition, including “meiotic recombination,” “RNA pol I transcription,” “peroxisome,” “oxidative phosphorylation,” “base excision repair,” and “RNA degradation” (Fig. 2B).

We recently developed BruUV-seq to map active enhancer elements based on UV-mediated repression of eRNA degradation by the RNA exosome.²⁶ Here we used BruUV-seq to assess the modulation of enhancer activity immediately following serum stimulation. Our results revealed that immediate changes in enhancer activity were common in neighborhoods of induced genes but not around suppressed genes. Of the 50 most highly induced genes, 18 genes were associated with 4 or more activated enhancer elements, 20 genes were associated with 1-3 activated enhancers, and only 12 genes were not associated with any activated enhancer element. In contrast, only one out of the 25 most highly downregulated genes was associated with an enhancer element with reduced activity. Thus, transcriptional induction immediately following serum stimulation is likely a result of enhancer activation, while immediate transcriptional suppression is likely achieved in an enhancer activity-independent manner, however further investigation of enhancer transcription and identification of enhancer target genes is necessary. It is likely that signal transduction pathways that respond to serum stimulation result in the activation of a set of transcription factors that primarily drive transcription initiation by binding to nearby enhancer elements. Transcription factor motif enrichment analysis indicated that SRF, CREB and ATF transcription factors may be driving the initial wave of induction (Fig. 3A). Interestingly, the top 3 most enriched motifs within suppressed genes did not have any identified transcription factor associated with them. This may reflect regulation of this understudied gene set by novel transcription factors.

If gene length plays an important role in coordinating expression following activation of particular pathways, we would predict that these gene length relationships are conserved across different species. Although genes tend to be larger in humans compared to other eukaryotes due to intron lengthening^7,31,33 we observed a strong correlation between gene lengths of orthologous genes in human and mouse. Therefore, we hypothesize that relative gene sizes may be important for the temporal regulation of gene expression in mammals. Though gene sizes may change over evolutionary time, the maintenance of relative gene sizes across species allows for a group of simultaneously induced genes to be expressed in a predetermined and presumably optimal order.

While we speculate that gene length acts as a biological timer for proper expression during the serum response (Fig. 9), the consequences of altering gene lengths and therefore the order of protein expression is unknown. The temporal order of gene expression during the serum response may only be important for certain sets of genes or for those involved in specific pathways. The immediately induced transcription factor genes showed a broad distribution of gene sizes, with gene lengths uniformly spread across a wide range as opposed to concentrated around a median value. Because this broad distribution may set up a distinct order of gene expression, we would predict that shuffling expression order would have significant biological effects. The fact that housekeeping genes are significantly shorter on average compared to all genes or transcription factor genes strongly supports the notion that the expansion of gene size of non-housekeeping genes has been selected for during evolution. Alternatively, introns in housekeeping genes have been shortened during evolution since these genes require constitutive rather than temporal expression. The larger sizes of serum-induced compared to serum-repressed genes implies that it is more important for cells to have a mechanism for temporally separating expression of genes that are induced compared to genes that are repressed. In addition to regulating the workflow at ribosomes, the temporal production of proteins as a consequence of different gene sizes may orchestrate their ordered assembly into multi-protein structures. Evolutionary selection of larger introns may also have been promoted by the advantage of placing regulatory elements and functional ncRNAs in the non-coding introns. Thus, introns can act as molecular timers with important biological regulatory functions both for protein production and assembly, and their size expansion has served as a critical element for evolutionary optimization.

In conclusion, our study provides a comprehensive profile of the early transcriptional response following serum stimulation. Early serum response genes display a wide range of lengths, which influences their transcriptional timing and allows for coordinated expression of functionally related genes. We speculate that gene size likely plays an important role in gene expression timing for many different cellular responses.

Materials and methods

Cell culture and serum stimulation

HF1, hTERT immortalized foreskin-derived human fibroblasts,^2,9,23-26,34 were grown in MEM supplemented with 10% FBS, L-glutamine, vitamin mix, and antibiotics. Starved cells were grown in MEM without the supplementation of FBS for 48 hrs. For serum stimulation experiments, FBS was added to the media of starved cells together with bromouridine to obtain 10% and 2 mM final concentrations of serum and bromouridine, respectively.

Bru-seq and BruUV-seq analyses

Bru-seq was performed as previously described.^24,25 Briefly, bromouridine (Bru) (Aldrich) was added to the media of starved or serum stimulated cells to a final concentration of 2 mM and incubated at 37°C for 30 min. BruUV-seq was performed as previously described ²⁶ using a UV dose of 100 J/m² prior to Bru-labeling. Total RNA was isolated using TRIzol reagent (Invitrogen) and Bru-labeled RNA was isolated from total RNA by incubation with anti-BrdU antibodies (BD Biosciences) conjugated to magnetic Dynabeads (Invitrogen) under gentle agitation at room temperature for 1h. cDNA libraries were generated from the isolated Bru-labeled RNA using the Illumina TruSeq kit, modified as described previously to produce strand-specific and size-selected libraries which were then sequenced using the Illumina HiSeq 2500 platform (50 base, single-ended reads).^24,25

Gene expression and serum response analysis

RPKM (reads per kilobase per million mapped reads) values were calculated for individual genes that were at least 300 bp long. For genes of lengths of 30 kb and less, RPKM values were calculated using read counts from the entire gene. For genes longer than 30 kb, an RPKM value was calculated using read counts from the first 30 kb downstream of the TSS. The R package DESeq³⁵ was used to test differential expression of genes whose mean RPKM between samples was greater than 0.5. Serum-response genes were defined as having a significant change in transcription initiation (adjusted p-values < 0.05, n = 2). Repressed genes were required to be expressed (RPKM value > 0.5) in the starved sample.

Ortholog size comparison

Mouse and human ortholog sizes were obtained using Ensembl Biomart.¹ The coefficient of determination (r²) was calculated. Genes with one-to-one orthologs were assessed.

Gene enrichment analysis

DAVID was used to analyze functional annotation enrichment.³⁶ We grouped induced genes into 5 categories based on gene length: 0-10 kb, 10-25 kb, 25-50 kb, 50-100 kb, and >100 kb. We grouped repressed genes into 2 categories: 0-15 kb and >15 kb. For the background gene list in each analysis, we used all expressed genes (RPKM > 0.5) which fell within the given size range. The pathways displayed have a p-value < 0.05. We also performed gene set enrichment analysis (GSEA).³⁷ Expressed genes were rank-ordered according to log2-fold changes compared to the gene expression in serum-starved cells. The top 10 pathways with the highest NESs were presented in Fig. 2.

Abbreviations

Bru	=	Bromouridine
CNV	=	Copy Number Variation
DAVID	=	Database for Annotation, Visualization and Integrated Discovery
GSEA	=	Gene Set Enrichment Analysis
hTERT	=	human Telomerase Reverse Transcriptase
MAPK	=	Mitogen-Activated Protein Kinase
qRT-PCR	=	quantitative Reverse Transcription Polymerase Chain Reaction
RPKM	=	Reads Per Kilobase of transcript per Million mapped reads
seq	=	sequencing
ST	=	Starved condition
SE	=	Serum-stimulated condition
UV	=	Ultra Violet light

Disclosure of potential conflict of interest

The authors have no conflict of interest with this study.

Acknowledgments

We thank the personnel at University of Michigan Sequencing Core for professional technical assistance and Manhong Dai and Fan Meng for administration and maintenance of the University of Michigan Molecular and Behavioral Neuroscience Institute (MBNI) computing cluster. We thank Dr. Artur Veloso and Dr. Thomas E. Wilson for development of bioinformatics analysis and visualization tools.

Funding

This work was supported by research grants from the National Institute of Environmental Health Sciences (1R21ES020946), National Human Genome Research Institute (1R01HG006786) and the National Institute of Health (P50CA130810). KSK was supported, in part, by an NIH training grant (T32GM07544) and BM was supported by the University of Michigan School of Public Health Environmental Toxicology and Epidemiology Program, National Institute of Environmental Health Sciences (T32ES007062).