3’UTR Shortening Potentiates MicroRNA-Based Repression of Pro-differentiation Genes in Proliferating Human Cells

In order for an APA event to potentiate targeting by miRNAs, potentially functional sites should exist 5' to the APA site, at a distance from the APA that is comparable to that typically seen between conserved miRNA sites and the canonical 3’ UTR ends. We have previously shown this region to be ~250 nucleotides from either end of the 3'UTR [12].

To identify transcriptome-wide APA sites, we subjected WI-38 human embryonic lung fibroblasts and their immortalized derivatives obtained through sequential serial transfers towards increased proliferation and transformation [13,14] to 3' sequencing and analysis [15]. The analysis included primary fibroblasts (“Control”), slow growers (early passage after immortalization), fast growers (extensive passaging after immortalization) and fast growers transformed by constitutively activated mutant H-RasV12 (“Ras”). 5765 genes were found to have at least one APA site in at least one of the cell types in this experimental system. Somewhat unexpectedly, we did not detect significant differences in the overall extent of global shortening between the different cell types. Remarkably though, when we aligned all genes with at least one APA event according to their most proximal APA site, a significant enrichment of conserved binding sites was observed within the 300 bases immediately upstream to the APA site (Fig 1B for control cells, S1A Fig for all time points). A similar picture emerged from the analysis of previously published 3'seq data [15,16], obtained in different cell lines (Fig 1C–1E). Of note, in those earlier experiments, APA was compared in MCF10A and BJ cells under growth arrest vs. proliferation and transformation, revealing significant shortening in the proliferating cells. Our inability to detect comparable shortening on a genome-wide scale in the WI-38 analysis may therefore stem from the fact that all compared states involved proliferating cells; alternatively, this may reflect the embryonic origin of WI-38 cells. Notwithstanding, though we do not see a genome wide tendency for shortening in these samples, we do observe many genes that do feature a shortened version, and we analyze them below.

A trivial reason for the similar enrichment of conserved miRNA sites upstream to APA sites in the different datasets could be that in all datasets the same genes undergo the same APA events. However, this pattern was still retained also when we performed a similar analysis only on APA-positive genes that differ between pairs of cell lines (examples in S1B–S1E Fig). Moreover, analysis of 3'seq data from mouse muscle tissue [16] revealed a similar peak of conserved miRNA binding sites upstream to the APA sites (Fig 1F). In the mouse tissue fewer genes were found to undergo APA, an observation which might reflect the highly differentiated state of the cells. Hence, in multiple cell types and in two mammals, distinct APA sites reside preferentially closely downstream to conserved miRNA binding sites, effectively repositioning such centrally-located miRNA binding sites and placing them in proximity to the 3’ end of the shortened transcript. In support of the emerging notion, when we compared the distribution of conserved miRNA binding sites between genes possessing at least one APA site in WI-38 cells and those without APA sites, we found (S2 Fig) that genes with APA sites tend to harbor fewer conserved binding sites near the distal end of their full length 3'UTR, relative to those without an APA (p-value = 6.5e-279, Student T-test). Conversely, genes with APA sites are relatively enriched in conserved miRNA binding sites within the proximal half of their 3’UTR (S2 Fig). This observation is intriguing: if APA merely serves to eliminate miRNA binding sites residing near the 3' end of the full length mRNA, one would expect APA-positive genes to be more enriched for functional miRNA binding sites near that end, providing them with an efficient on/off switch controlled by APA. The fact that the opposite trend is actually observed strongly suggests that the interplay between APA and miRNAs may allow regulation that is richer than mere binding sites elimination. Specifically, this may serve as further indication that the shorter 3’ end, positioned immediately upstream to the APA site, can dynamically potentiate new miRNA binding sites as they become positioned closer to the 3’ end of the shorter transcript.

Analysis of the conservation state of each miRNA binding site by itself is an accepted indication for its functionality [10,17,18]. However, such conservation might be due to other attributes, for example another functional feature of the 3’UTR residing in this location. We therefore asked whether the high conservation of the miRNA sites located upstream to APA sites can be attributed to conservation of their neighborhood, or is preferentially targeting the miRNA sites? To address this question we looked at the conservation profile around APA sites, for genes with and without conserved miRNA binding sites in the 300 bases 5’ to the APA sites. As can be seen in S3 Fig, the profile of the genes without conserved sites 5’ to APA sites is lower specifically in this area, indicating that the high conservation present for the other genes might be due to the presence of miRNA targeting, and not merely APA sites. We also looked at the conservation profile of all conserved miRNA binding sites located in the 300 bases 5’ to APA sites. For that we used PhastCons and PyhloP [19,20] as provided by the UCSC browser. We observed a sharp peak of conservation exactly overlapping the 7 bases of the binding sites, whereas its surroundings are significantly less conserved (Fig 2A), strongly arguing in favor of a selective pressure to conserve specifically the conserved miRNA binding sites. Another measurement of site conservation is the P_CT score, which controls for the 3’UTR surroundings, dinucleotide conservation and other parameters unrelated to miRNA functionality [21]. Importantly, this score allows to assess the extent to which conservation of a site is likely to be due to miRNA functionality. We compared the distribution of P_CT scores for miRNA binding sites 300 bases 5’ and 3’ to APA sites, and found that the scores are significantly higher for the sites located just 5’ to the APA sits (Fig 2B, p-value = 2e-7, Student T-test). Moreover, we employed the Context++ scoring system of miRNA binding sites, which takes into account many additional parameters of each site and its surroundings beyond mere conservation, and provides a score for the probability that this site is indeed functional [22]. We compared the scores of the binding sites 300 bases 5’ and 3’ to APA sites, for conserved and non-conserved sites. For both groups, the scores for sites located 5’ to APA sites were significantly lower (hence indicating higher functionality) (Fig 2C and 2D, p-value = 2e-6, 7e-165, Student T-test). Together with the P_CT score analysis, this argues that these sites are more likely to be indeed functional.

In sum, the above analyses strongly suggest that the miRNA binding sites located just 5’ to APA sites are functional beyond sequence conservation. Since the Context++ scores of non-conserved miRNA binding sites were also better for the ones that are positioned just 5’ to APA sites, we can conclude that even the non-conserved binding sites are probably more functional for miRNA targeting when located in that region.

In non-proliferating cells, most genes express mainly the full length version of their 3’UTR. However, in proliferating cells and particularly in cancer, increased usage of APA sites has been shown to enable miRNA binding sites elimination from the 3’UTRs of proliferation-associated genes [1,6]. By the same rationale, we predicted that if our above observations are physiologically relevant, then genes becoming more susceptible to miRNA regulation (and hence more effectively repressed) owing to APA would tend to be those that should be preferentially downregulated during proliferation and in cancer. To address this prediction, we initially compared two gene sets: those at the core of the cell cycle machinery, and those involved in patterning of the embryo during development. We chose these two gene sets as they represent two opposing classes of archetypical proliferation and differentiation genes [23]. Yet, these gene sets are relatively small and for only a portion of them we detected APA events. We therefore sought to expand these gene sets to include functionally related genes and thus gain further statistical power. To that end, we expanded each gene set to include additional genes either by similar codon usage (see M&M) or based on the GSEA tool [24]. We then analyzed the miRNA binding site landscape around the APA sites of the two gene sets. Remarkably, we observed a significant difference between the two groups of genes: while the proliferation-related genes (“Pro-Prolif.”) display only modest enrichment of conserved miRNA binding sites immediately upstream to the APA site, the differentiation-related genes (“Pro-Diff.”) show a markedly elevated abundance of conserved miRNA binding sites in the corresponding region (Fig 3A and S4A–S4C Fig). This strongly suggests that differentiation-related genes are more prone than proliferation-related genes to regulation by miRNA binding site potentiation via APA. This is in line with the documented increased APA usage during proliferation and cancer, when differentiation-related genes are expected to be downregulated.

These results suggest that, in addition to its documented ability to alleviate miRNA-mediated repression of proliferation genes, 3’ UTR shortening is also used to potentiate preferentially the repression of differentiation genes.

To obtain clues about the miRNAs that bind sites potentiated by APA, we performed comparative miRNA microarray analysis on WI-38 cells and their progressively transformed derivatives [13,14]. We found that most of the miRNAs upregulated particularly in the highly proliferative stages (fast growers and Ras-transformed) belong to the miR-17-92 cluster (Fig 3B and S5 Fig). Indeed miR-17-92 is a well-studied proliferation-associated cluster [25], which includes 6 miRNAs with 4 different binding site sequences. Target genes of these miRNAs are significantly enriched for regulation of differentiation and negative regulation of proliferation (S1 Table), consistent with the notion that these miRNAs promote cell proliferation. Notably, in comparison to all miRNAs on the array, conserved sites for members of the miR-17-92 cluster are significantly enriched immediately 5' to APA sites (Fig 3C). Furthermore, while for all miRNAs in the genome we see a marked increase in the abundance of conserved binding sites near the 3' end of the full length 3’UTR, conserved binding sites of the miR-17-92 cluster are relatively less enriched in that region (Fig 3D). Thus, as compared to the bulk of the cellular miRNAs, miR-17-92 cluster members preferentially bind sites that are located upstream to APA sites rather than near the distal end of the full length 3’UTR. APA is therefore expected to preferentially potentiate the repressive effects of those proliferation-associated miRNAs.

Overall, the above findings further support the conjecture that while APA enables proliferation-associated genes to escape miRNA regulation, it confers increased regulation upon pro-differentiation genes.

Low passage WI-38 cells were obtained from ATCC (CCL-75). Slow growing hTERT-immortalized WI-38 cells and their fast growing derivatives, generated by extended passaging in culture of the slow growing cells [13], were kindly provided by Varda Rotter. Ras-Transformed derivatives of the fast growing WI-38 cells were obtained by infection with a recombinant retrovirus expressing H-RasV12 as described in [14].

WI-38 cells were grown in 37° in MEM supplemented with 10% non-heat-inactivated fetal bovine serum (Sigma), pen-strep, sodium pyruvate, L-glutamine solution (Beit HaEmek). RNA was extracted using Nucleospin miRNA kit (Macherey-Nagel), according to manufacturer’s instructions.

miRNA array analysis was done in duplicates, using the miRNA Complete Labeling and Hyb Kit (Agilent, 5190–0456) according to the manufacturer's instructions. Briefly, for each sample 100ng RNA was dephosphorylated, denatured, labeled with Cyanine 3-pCp and purified using Micro Bio-Spin 6 Columns. Hybridization was done for 20h with Agilent SurePrint G3 Unrestricted miRNA 8x60K (Release 19.0) arrays. Arrays were scanned using an Agilent DNA microarray scanner, and analyzed using the AgiMicroRna package in R [26] with the RMA algorithm. The heat map was generated with Matlab Clustergram function for the mean fold change of each sample vs. control (mean of the duplicate arrays).

The 3’Seq protocol is based on Jenal et al. [16] and incorporates additional modifications described in Martin et al., 2012[27]. Basically, 25 μg of total RNA were heat-fragmented for 12 minutes in 1x Fragmentation Buffer (Ambion) at 70°C to generate RNA fragments of ~100 nucleotides. Next, the 3’end poly(A) RNA fragments were selected using the Oligotex mRNA Kit (QIAGEN) and RNA was end-repair with T4 polynucleotide kinase for 45 minutes at 37°C following manufacturer’s instructions. Afterwards, RNA 3’ends were blocked for ligation by incubation with 1mM Cordycepin 5′-triphosphate (Sigma) and 10U of polyA polymerase (PAP, NEB), in 1x PAP buffer for 30 minutes at 37°C. Finally, a P7 RNA adapter (5’-CAAGCAGAAGACGGCAUACGAGAU-3’) was ligated to the 5’end using 2U of T4 RNA ligase I and 2.5uM of RNA adapter, for 4h at room temperature. Between each step, RNA was purified using Agencourt RNAClean XP magnetic beads (Beckman Coulter) following the manufacturer’s instructions. At this point, RNA fragments were converted to cDNA employing the Superscript III RT kit (Life Technologies) and an anchored oligo(d)T stem loop primer containing a barcoded Illumina adaptor as in Martin et al (See S3 Table). Next, cDNA was purified twice with Agentcourt AMPure XP magnetic beads (Beckman Coulter) using a ratio 1.5:1 beads:sample. To generate the final 3’Seq library, the cDNA with the correct adaptor sequences was enriched/amplified using Phusion DNA polymerase (Life Technologies) and primers P7 and Illumina_Truseq, for 12 cycles following manufacture’s recommendations. Finally, the 3’seq library was size selected with AMPure XP magnetic beads by two rounds of purification with a ratio 1:1 beads:sample, before being sequenced on an Illumina HiSeq2000 system.

In our protocol of profiling transcript 3' ends, sequenced reads start with a barcode for sample multiplexing which is followed by six Ts whose end marks the precise location where the poly(A) tail starts. These six Ts therefore allow the mapping of poly(A) cleavage sites (CSs) with a nucleotide resolution. After trimming the barcode and six Ts, reads were aligned to the human genome (hg19) using TopHat [28]. Up to two mismatches were allowed in the reads’ seed region (the first 28 nt). As CS location often fluctuates around a major site, we merged reads from all samples and identified "read runs” (that is, genomic intervals that are "tiled" by multiple reads where the distance between the start of consecutive ones is below 10 nt). We considered the local maxima of these runs as the CS locations. We required a spacing of at least 50 nt between consecutive CSs. (In case of lower spacing between CSs, the one supported by a higher number of reads was chosen). Only CSs supported by at least 10 reads (at the location of the CS run maximum) were considered in subsequent analyses. The median length of the runs was 11 nt. Overall, 41,972 CSs were detected in our dataset. Priming of the oligo-dT primer to genomic regions that are A-rich (“internal priming”) could lead to false call of CSs. To reduce the rate of such false calls we extracted genomic sequences of 50 nt centered at the location of the putative CSs, and filtered out CSs that contained in that region a stretch of 10 nt of which at least 8 were As and the rest were Gs. 4,307 suspected CSs were filtered out.

miRNA binding sites were defined as perfect 7-mers, which are reverse complement to the seed of the miRNAs, for all human and mouse miRNAs listed in miRBase release 17 [29]. Conserved binding sites were taken from TargetScan release 6.2 [18].

P_CT scores of all miRNA binding sites of type “7mer-m8” (perfect 7mer) for the conserved miRNA families were taken from TargetScan release 6.2 [18]. Context++ scores of all miRNA binding sites of type “7mer-m8” were taken from TargetScan release 7 [22]. Only miRNA binding sites from genes with at least 500 bases around the APA sites were taken into the analysis.

All analyses were done only for genes whose accession number was included in TargetScan 3’UTRs list (as defined in their website). A cleavage site was assigned to a gene if it was included in the coordinates of its 3’UTR (and 20 bases further, after the 3’ end of the 3’UTR). A gene was considered to have an APA site only if it has at least two cleavage sites assigned to its 3’UTR. In all analyses, the APA site that we took into account was the 5’ most in the 3’UTR.

To compute the statistical significance of the main signal—enrichment of conserved miRNA binding sites at particular distances from the APA site, we created a randomized null model. In this null model each binding site’s location was recorded and for each gene a random position of an APA site was drawn from the full length 3’ UTR (omitting the first and last 1000 bases in order to allow inspection of that vicinity around the randomly chosen location). The distance from the binding site and the randomized APA location was computed and the procedure was repeated 10,000 times. This yielded a distribution of distances as shown in all plots, once for conserved binding sites, and once for non-conserved binding sites. * indicates p-value < 0.05 for the null hypothesis that for a specific distance from the APA site, the number of conserved binding sites is similar as in a random APA site, or higher.

PhastCons and PhyloP scores of each base in the genome (hg19) were taken from UCSC [19,20]. The profile around APA sites of genes was for genes with at least 1000 bases from each side of the APA sites. The profile around conserved miRNA binding sites was for all miRNA binding sites located in the 300 bases 5’ to the APA site. The sites were aligned and the mean conservation profile was calculated.

We defined the “Pro-Proliferation” and “Pro-Differentiation” gene sets as follows: We began with the Gene Ontology sets termed “M-phase of cell cycle” and “Pattern Specification”, two gene sets that were recently shown [23] to serve as archetypical proliferation and differentiation genes, with distinct codon usage. To augment the number of genes belonging to each of the two sets we searched the entire genome for additional genes whose codon usage was highly similar to either of the two groups, thus expending the two sets from 92 and 82 genes originally to 229 and 136 (For correlation threshold 0.7: 226 and 147, for correlation threshold 0.8: 172 and 125). We expanded the two original gene sets (“M-phase of cell cycle” and “Pattern Specification”) also by correlation to other GO groups. These groups were taken as the top 50 groups for each original group using the GSEA tool track C5 [24].

After computing the miRNA binding site distribution around APA sites for the genes in each of the sets we estimated a p-value on the difference between conserved miRNA binding sites density at each distance from APA sites as follows: we repeated 10,000 on randomly partitioning the genes with high correlation to the two groups into two groups, one with 229 genes and one with 136 genes for the codon usage expansion threshold 0.75 correlation, and 175 and 167 genes in the expansion by the GSEA tool. In each such random partition we recorded, at each location relative to the APA, the fraction of conserved miRNA binding sites. The p-value was estimated as the fraction out of the 10,000 repetitions in which the real partition into “Pro-Proliferation” and “Pro-Differentiation” resulted in a difference in binding sites count. * indicates one-sided p-value < 0.05 for the null hypothesis that for two groups of these sizes, the difference in number of conserved binding sites is as for the two original groups or higher.

The miR-17-92 binding sites analysis was similar to the codon usage. Here too randomization was done 10,000 taking a random group of genes in the same size of the genes with binding sites for miR-17-92 miRNAs (104), and computing for each distance from the APA site the difference in the number of conserved binding sites for the random group and for all genes. * indicates one-sided p-value < 0.05 for the null hypothesis that for a random group of genes in the same size as the original one, the difference between the number of conserved binding sites between this group and all genes is as good as for the original group or higher.

Analysis of target genes for miR-17-92 cluster miRNAs was done using Gorilla tool [30]. Genes with at least two conserved binding sites for miRNAs in the cluster were defined as the “target set”, and were analyzed in comparison to the background group which was the whole genome.

All relevant data are within the paper and its Supporting Information except for microRNA array files, which are available from: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE77338↗.