A quantitative atlas of polyadenylation in five mammals

We developed PolyA-seq, a strand-specific and quantitative method for high-throughput sequencing of 3′ ends of polyadenylated transcripts, and used it to globally map polyadenylation (polyA) sites in 24 matched tissues in human, rhesus, dog, mouse, and rat. We show that PolyA-seq is as accurate as existing RNA sequencing (RNA-seq) approaches for digital gene expression (DGE), enabling simultaneous mapping of polyA sites and quantitative measurement of their usage. In human, we confirmed 158,533 known sites and discovered 280,857 novel sites (FDR < 2.5%). On average 10% of novel human sites were also detected in matched tissues in other species. Most novel sites represent uncharacterized alternative polyA events and extensions of known transcripts in human and mouse, but primarily delineate novel transcripts in the other three species. A total of 69.1% of known human genes that we detected have multiple polyA sites in their 3′UTRs, with 49.3% having three or more. We also detected polyadenylation of noncoding and antisense transcripts, including constitutive and tissue-specific primary microRNAs. The canonical polyA signal was strongly enriched and positionally conserved in all species. In general, usage of polyA sites is more similar within the same tissues across different species than within a species. These quantitative maps of polyA usage in evolutionarily and functionally related samples constitute a resource for understanding the regulatory mechanisms underlying alternative polyadenylation.

Sequencing of mRNA and noncoding RNA has made important contributions to our understanding of biology and disease, with numerous implications for diagnostics and therapeutics. As an outcome of rapidly expanding sequencing capabilities, recently described methods produce comprehensive representations of the transcriptome (Mortazavi et al. 2008;Armour et al. 2009;Wang et al. 2009;Levin et al. 2010) and have been used to discover and monitor alternative splicing (Sultan et al. 2008;Wang et al. 2008;Wilhelm et al. 2008), as well as gene expression (Marioni et al. 2008) and its underlying regulatory genetic variation (Montgomery et al. 2010;Pickrell et al. 2010). While transcriptome sequencing studies continue to focus on gene expression and RNA processing, mapping of polyA sites has received considerably less attention, despite evidence suggesting that alternative polyadenylation is common in metazoans (Lee et al. 2007;Ozsolak et al. 2010) and contributes to phenotypic variation and disease. Avoidance of microRNA regulation via alternative polyA sites, for example, plays a role in development (Mangone et al. 2010; Thomsen et al. 2010;Jan et al. 2011) and cancer Mayr and Bartel 2009). Furthermore, extensive usage of tissue-specific sites, some of which are associated with cis-regulatory elements, suggests that alternative polyadenylation is tightly regulated (Proudfoot et al. 2002) and has important physiological implications. Lastly, the 39UTRs of some genes are expressed independently of the cognate coding regions, and may thus possess previously unappreciated regulatory functions (Mercer et al. 2011).
A simple and reliable method for high-throughput mapping of polyA sites is needed to help define the complete set of transcripts present in a dynamic transcriptome and to better understand the role of alternative polyadenylation. ESTs and full-length cDNAs have traditionally been used for defining primary transcripts, but become expensive and laborious on a large scale; consequently, the genomes of many species remain sparsely annotated. Highthroughput variations on serial analysis of gene expression (SAGE) (Velculescu et al. 1995) have proven useful in DGE applications (Asmann et al. 2009), but their reliance on restriction sites precludes precise mapping and excludes some genes entirely. Paired-End diTags (PETs) simultaneously capture 59 and 39 ends of mRNAs (Ng et al. 2005(Ng et al. , 2006Fullwood et al. 2009), although the short tag sequences and complex cloning protocol impose practical challenges for genomic alignment and assay scalability.
Recently, several studies have described methods that leverage high-throughput sequencing to map 39 transcript termini. The earliest was a modified RNA-seq protocol, where reverse transcription is primed with a T 20 VN oligo (20 Ts followed by a non-T, then a random base) and the second strand is made with RNase H + DNA Pol I, followed by standard library construction (Yoon and Brem 2010). Using paired-end sequencing, one mate captures sequence near the end of the transcript, while a stretch of either As or Ts in the other mate reveals the orientation relative to genomic DNA. In a method called MAPS, Fox-Walsh et al. (2011) use a biotinylated variant of T 20 that enables bypassing mRNA isolation from total RNA and produces multiplexing-enabled 39-biased cDNA libraries. The main limitations of both of these approaches are that the majority of sequencing reads do not precisely map to the transcript termini and libraries are contaminated with internally primed fragments. Jan et al. (2011) developed an elegant molecular approach that minimizes representation of internal priming events by requiring ligation of a splint-oligo to the end of the polyA tail. The method was used to reveal thousands of novel polyA sites in C. elegans, but requires many steps during library construction and has not been shown to be quantitative. Lastly, PAS-seq is similar to MAPS in that a universal sequence is upstream of the T 20 VN, which enables direct amplification by serving as a PCR anchor (Shepard et al. 2011), but enables precise definition of the polyA site by sequencing directly into the end of the transcript. All approaches were thus far applied to small numbers of samples.
We developed a novel method that is similar to PAS-seq (Shepard et al. 2011), but differs in its use of random priming for second-strand synthesis, enabling rapid library generation. Our protocol takes 6 h, of which only 2 h are hands-on. Instead of T 20 VN, we used a T 10 VN, which also facilitates the use of a custom sequencing primer with more favorable hybridization kinetics during sequencing. Our method, called PolyA-seq, captures transcript sequence immediately upstream of polyA stretches, retains strand specificity, and accurately represents transcriptional abundances as demonstrated by agreement with MAQC data (Shi et al. 2006). Thus, it also quantitatively captures the usage of any given polyA site. We developed a robust filtering scheme to remove internal priming events and calibrated the false-discovery rate through experimental validation. Here, we describe the method and an exploratory analysis of a polyA site atlas that we generated in five mammals, focusing on tissue and species specificity to gain insight into general evolutionary trends in polyadenylation.

Results
PolyA-seq is a simple high-throughput strategy for sequencing 39 transcript termini Library construction consists of the following steps: (1) reversetranscription primed with an oligonucleotide consisting of a universal sequence that serves downstream as a PCR anchor, followed by T(10)VN (i.e., 10 thymidines, then a random base other than a thymidine, then a random base), (2) second-strand synthesis using random hexamers linked to a second PCR anchor, and (3) nested PCR to add Illumina-specific adapters while preserving strand orientation ( Fig. 1A; See Supplemental File 1 for detailed protocol). A primer ending in 10 Ts is then used for sequencing. This protocol yields an antisense read beginning with the base immediately upstream of the polyA stretch. Amplicons are typically 200-500 bp after accounting for adapter sequences (Fig. 1B) and are compatible with standard cluster generation protocols. We applied PolyA-seq to 24 samples from five species, with each sample sequenced in one lane on an Illumina GAIIx sequencer. We aligned reads to the genome, transcripts, and exon-exon splice junctions, and retained all uniquely mapped reads for further analysis ( Fig. 1C; see Methods for more details).

PolyA-seq is an accurate approach for DGE
To assess the quantitative potential of PolyA-seq, we assayed the Human Brain Reference and Universal Human Reference (UHR) MAQC samples (Shi et al. 2006) in technical replicates. We measured expression as reads per million, namely, the number of reads aligned to a transcript divided by the total number of uniquely aligned reads in the experiment (in millions). Since we were interested in demonstrating that PolyA-seq is a simple approach for DGE applications, we used all mapped reads for this analysis, rather than filtered reads (see below); however, using filtered reads had a negligible impact on DGE accuracy (data not shown). We found reproducibility to be high (Brain technical replicate Pearson r = 0.994, Fig. 2A; UHR technical replicate Pearson r = 0.988). We also observed good agreement with qRT-PCR (Pearson r = 0.948 for Brain/UHR ratio; Fig. 2B). In a systematic comparison of six expression profiling approaches applied to MAQC data (qRT-PCR [Shi et al. 2006], Affymetrix microarrays [Shi et al. 2006], Agilent microarrays [Shi et al. 2006], RNA-seq [Bullard et al. 2010], and 39-tag DGE [Asmann et al. 2009]), PolyA-seq performed at least as well as any other method in both ratio and absolute gene-expression quantification ( Fig. 2C; see Methods for details). RNA-seq and PolyA-seq attained the highest agreement with qRT-PCR in the ratio-based comparisons. Motivated by the possibility that one of these approaches may be superior when sequencing throughput is limited (e.g., when multiplexing), we assessed correlation with qRT-PCR with expression values computed from randomly selected subsets of mapped reads. Both methods were considerably less accurate with fewer than 1 million mapped reads, but did not differ significantly in their performance at any level of input (Fig. 2D). Both approaches yielded similar proportions of uniquely mappable reads (71% for RNA-seq; see Table 1 for PolyA-seq).

Discriminating genuine polyA sites from internal priming events
Since PolyA-seq relies on priming from a stretch of adenines, it captures polyA tails added post-transcriptionally, as well as internal (i.e., transcribed) stretches of adenines encoded in the genome. EST, full-length cDNA, GIS-PET (Fullwood et al. 2009), PAS-seq (Shepard et al. 2011), and polyA capture cloning and sequencing (Mangone Input was polyA+ selected RNA (green). Reverse transcription using U1-T10VN was followed by RNase H treatment to degrade RNA. Second-strand synthesis using U2-N6 was achieved through a random-primed Klenow extension. U1 and U2 have sequence complementarity to Illumina-specific adapters, which are added through PCR. This yields DNA libraries that can be directly sequenced. (B) A typical library consists of amplicons ranging from 200 to 500 bp (Illumina adapters account for 79 bp). (NTC) No-template control. (C ) Computational procedure: reads were aligned to the genome and transcriptome ([*] defined here as known and predicted splice junctions extracted from UCSC Known Genes, RefSeq, and Ensembl, followed by conversion to genomic coordinates; see Methods for more details). Matches with unique loci were then filtered on internal priming potential and clustered into polyA sites. et al. 2010) strategies also exhibit some degree of internal priming. Building on previous efforts (Lee et al. 2007), we used genomic sequence immediately downstream from our reads to distinguish genuine polyA sites from internal priming events. We built an empirical model that estimates the probability that any given site reflects a genuine polyadenylation event based on the 10 bases immediately downstream from the site (corresponding to the length of the RT primer). To build the model, we first compiled lists of internal priming events and genuine polyA sites experimentally. We generated PolyA-seq libraries from UHR, but used an oligo ending with 10 Ts (i.e., no VN) for first-strand synthesis. This led to random placement of PolyA-seq reads within A-stretches. A total of 923,911 reads contained genome-alignable 59 ends, but ended with unalignable stretches of As, thereby defining 28,818 unique genuine polyA sites (see Methods for further details). Reads whose 39 ends could be perfectly aligned to the genome were treated as internal priming events (see Methods and Supplemental Figure 1 for more details). Some of these reads may, in fact, represent polyA sites, since the T10 primer can hybridize to the 59 end of a polyA tail by chance. However, we estimate this proportion to be small, and these false-negatives can only lead to an underestimate of false discovery (see Methods).
We found a striking difference between priming site base frequencies at internal and polyA sites (Supplemental Fig. 1), and used these distributions to construct a log-odds scoring model that estimates the probability of a site being genuine vs. internally primed. Leveraging this model and the map of internal priming and polyA events in UHR, we calibrated our filtering threshold to a sensitivity of 85.6%. In the UHR data set this corresponds to a false-discovery rate (FDR) of 2.5% (Supplemental Fig. 2; see Methods) that was also consistent with additional RT-PCR based validation . Given the similar genomic and transcriptional complexities of all samples sequenced (Chan et al. 2009), we expect this number to be generally representative of the FDR in mammalian samples. Base frequencies of filtered sites strongly resemble those of well-established gene termini (Supplemental Fig. 6A), and our model is substantially better than simple filtering based on the number of downstream adenines as implemented previously (Lee et al. 2007;Shepard et al. 2011) (Supplemental Table 1; e.g., Supplemental Fig. 6B). Furthermore, we observed hundreds of previously reported polyA sites that are likely internal priming events (e.g., Supplemental Fig. 6C). We note that simply using a sequencing primer that lacks 10 terminal Ts to identify genuine polyA sites would not reveal the original RNA priming site, since the primer used in firststrand synthesis is incorporated into the sequencing library (not the RNA priming site) (Fig. 1A). Finally, we note that including a parameter that represents the read counts for each site to aid in discriminating real sites from internally primed sites did not improve performance (data not shown). We reasoned that real sites may have overall higher usage (i.e., more reads supporting them), but this was not the case.
From visual inspection of the data it was apparent that known polyA sites typically corresponded to a major cluster of reads mapped to a single position, with several minor peaks within a few bases. This ''wiggle'' is also evident in EST and cDNA alignments (Fujita et al. 2011) and is a likely consequence of the inexact nature of transcriptional cleavage (Proudfoot et al. 2002). We therefore clustered all peaks within 30 bases on the same strand, and retained only the highest peak within each cluster; varying this window from 10 to 200 bp did not significantly impact the number of sites (Supplemental Fig. 7). Table 1 summarizes the number of PolyA-seq reads and polyA sites that we identified in each sample. Reproducibility of PolyA-seq is high; 83% and 84% of Brain and UHR polyA sites detected in the MAQC samples, respectively, were also identified in the corresponding technical replicate.
PolyA-seq is precise and captures many known, novel, and alternative polyA sites in all species We evaluated the precision and complexity of our atlas in human, where the transcriptome is extensively annotated (compared with rhesus, for example), and where we had the most data (9/24 samples). A total of 99.3% of UHR PolyA-seq reads that overlapped with a known RefSeq 39UTR (excluding 8% of 39UTRs that overlapped other genes on opposite strands) agreed on strand of origin, demonstrating the high strand specificity of PolyA-seq (example shown in Fig. 3A). The vast majority of our human polyA sites agree with known transcript termini to single-base precision ( Fig.  3A-C). The single dominant peak in Figure 3B corresponds to known transcript termini (see Methods for further details), and the sites situated upstream are likely alternative polyA sites (tissue-dependent example shown in Fig. 3C), which, in human, outnumber uncharacterized transcript extensions (Fig. 3B).
To begin to explore the distribution of alternative polyA sites, we first clustered 39UTRs annotated by RefSeq (Wheeler et al. 2003), UCSC Known Genes (KG) (Fujita et al. 2011), and Ensembl into 27,175 unique 39UTR models (see Methods for further details). We also extended each model by 1 kb to capture uncharacterized transcript extensions. We then asked how many sites are detected per UTR model in each human tissue. To enable comparative analysis we accounted for different numbers of aligned input reads by subsetting each human sample to 9,886,234 aligned reads (to match human muscle) by random selection. A total of 20,873 (76.8%) 39UTR models were detected in at least one human tissue; on average, 46.3%, 27.9%, and 30.7% of UTRs had one, two, and three or more detected sites (Fig. 3D). Aggregating reads across tissues revealed a shift toward more sites per gene, indicating that many of these sites are tissue-dependent (Fig. 3C). Overall, dominant (most utilized) sites account for >90% reads (Fig. 3D), suggesting that, although prevalent, secondary sites tend to be used much less frequently than dominant sites.
In all five species we observe excellent concordance with known polyA sites. Most PolyA-seq reads map to known polyA sites, although most polyA sites are novel, particularly in rat, dog, and rhesus, where existing annotations are not as extensive as in human and mouse ( Fig.  3G; see Methods for further details). This result is consistent with our observation regarding usage of dominant sites, suggesting that novel sites are often used at lower levels, and may explain why they have not been previously reported. It is unlikely that low-frequency sites are enriched for false-positives, since inclusion of a parameter that captures site usage in our filtering model did not improve performance. In other words, low-frequency sites are just as likely to be real sites as are high-frequency sites. A primate-specific polyA site in a microRNA cluster is shown in Figure 3F, occurring 75 bp upstream of a conserved site.

Polyadenylation in noncoding RNAs
In addition to detecting polyadenylation in mRNAs, PolyA-seq captures polyadenylated noncoding RNAs, including primary microRNA transcripts and antisense RNAs. For example, polyA sites are detected downstream from numerous let-7 and other microRNA clusters in all or nearly all tissues, corroborated by ESTs in human, mouse, and rat, but not dog or rhesus, where ESTs are scarce (e.g., Fig. 3F; Supplemental Fig. 8). Polyadenylation downstream from the liver-specific microRNA MIR122 is detected very highly in liver but not in any other tissue (Supplemental Fig. 9A); similarly, MIR124-1, known to be brain-specific, is primarily detected in brain (Supplemental Fig. 9B). In species with abundant EST data, such as human and mouse, ESTs corroborate PolyA-seq sites. In other species, particularly dog and rhesus, PolyA-seq supplements predictions based on sequence conservation by providing direct evidence of transcription. Noncoding RNAs detected via PolyA-seq also include the human accelerated region (HAR) transcripts HAR1A and HAR1B (Supplemental Fig. 10A), XIST (Supplemental Fig. 10B), and HOTAIR (Supplemental Fig. 10C), as well as transcripts antisense to DLX1 (Supplemental Fig. 11A) and HOXA11 (Supplemental Fig. 11B).

Motif enrichment analysis
In each tissue, and for all possible 6-mers, we assessed the frequency of 6-mer occurrence within 150 bp of all polyA sites and compared it with that observed in randomly selected 39UTR sequences of the same length ( Table 2). The top 100 most statistically enriched 6-mers were identified by calculating Z-scores for the difference between tissue and background using a normal approximation to the binomial distribution. Positional bias for each k-mer was then assessed by calculating a x 2 statistic assessing the null hypothesis of uniform positional distribution of the k-mer within 150 bp of polyA sites in a tissue. The top 10 positionally  . In human and rhesus, two polyA sites (purple spikes) correspond to two canonical polyA signals (AATAAA; black tick marks), the first of which is present only in primate genomes (data not shown). In rat, mouse, and dog, only the downstream polyA site is detected, in accordance with the absence of the upstream polyA signal. (G) Distribution of reads and polyA sites across genomic features. All reads were aggregated in each species and then filtered and clustered as described in the main text.
enriched k-mers from each tissue, as well as the frequency of occurrence of the 6-mer at any position near a polyA site are provided. In all cases, these frequencies are highly enriched relative to background (Fisher's exact test P-values below the level of computing precision in all cases). A total of 98% of all polyA sites identified contained a match to one of these hexamers and we observed the expected enrichment upstream of sites (Fig. 4A). Individual base frequencies (Fig. 4B) also supported the known CA dinucleotide enrichment at the cleavage site (Sheets et al. 1990;Chen et al. 1995).

Evolutionary conservation of polyA site usage
To our knowledge, this study is the first systematic and unbiased effort to compile polyA sites in matched tissues across multiple species in parallel. Given the extent of tissue-specific alternative polyA-site usage, an obvious evolutionary question is whether polyA-site usage is conserved across species or whether species carry unique polyadenylation signatures that may drive phenotypic variation. We reasoned that if polyA-site usage was in general under purifying selection, usage should be more conserved in the same tissue across species than in different tissues within the same species.
To begin to address this question we compiled a five-way orthologous set of genomic coordinates for each polyA site in our atlas. A total of 2590 polyA sites were detected in at least one tissue of every species (Fig. 5A). We normalized site usage in each sample to Z-scores (number of standard deviations away from the mean number of reads/site) and clustered sites ( Fig. 5B; hierarchical, Pearson correlation-based distance, average linkage). To comprehensively represent the similarities among samples, we also clustered the Pearson similarity values (Fig. 5C). We found that polyA-site usage is often strikingly similar in orthologous tissues, particularly in brain, kidney, and liver (Fig. 5C). Clusters of similar tissues are also better pronounced between closely related species. Rhesus and human tissues, for example, as well as mouse and rat tissues, tend to be more strongly correlated. Conservation of gene expression across similar tissues has been shown previously (Chan et al. 2009) and may be contributing to the general conservation of polyadenylation that we observe. In any case, this observation suggests that alternative polyadenylation is not a dominant force driving phenotypic diversity. Differences in conservation patterns among tissues present an opportunity for future research.

Discussion
Our primary intent was to generate an accurate, comprehensive, and systematic mapping resource to augment gene structure annotation efforts and to empower additional functional studies. The union of the top ten hexamers from five tissues are shown along with a ranking within each tissue. Following our robust filtering scheme, experimental validation revealed high sensitivity (;85%) and low false discovery (;2.5%), and strong agreement with previously defined polyA sites ( Fig. 3; Table 1). Positional enrichment for the canonical polyA signal (Fig.  3) supports these strong performance metrics. We discovered more than 60,000 novel polyA sites in each of the five species, and more than ;150,000 in rat, dog, and rhesus, where previous annotation is limited (Table 1). In all species the ratio of novel to known sites is greater in the cumulative tally of unique sites than in any individual tissue, suggesting that most of the novel sites are tissuespecific. While some genes appear to predominantly express a single polyA site across normal tissues and others harbor tissuedependent sites (e.g., ELAVL1; Supplemental Fig. 6B), close to 70% of genes consistently use multiple polyA sites across all tissues in which they are expressed.
Although the majority of sequencing reads support known polyA sites, most sites, particularly in species where annotation is extensive (e.g., human, mouse), are novel (Fig. 3G). This implies that polyadenylation at novel sites tend to be low-frequency events (i.e., supported by few reads). Our validation/filtering efforts strongly suggest that these sites are genuine polyadenylation events (see Results), although their physiological roles are unclear. Similar to early reports of widespread transcription (Cheng et al. 2005), it was later shown that by mass the transcriptional origin of most nonribosomal RNA in a cell can be traced to coding exons or pre-mRNAs (van Bakel et al. 2010). Indeed, the functional roles of thousands of novel transcripts were widely debated, although it was later shown that a substantial class of ''novel transcripts'' plays more general roles in transcriptional regulation (Neil et al. 2009;Xu et al. 2009). Similarly, we believe that there will be examples of novel low-frequency polyA sites that have important functional implications, but the majority are likely polyadenylation events of cryptic transcripts or aberrant events that are tolerated by the cell. Nonetheless, this atlas contains thousands of novel high-frequency sites that will aid transcriptional annotation efforts, particularly in less defined species (e.g., rhesus, dog).
PolyA-seq provides concrete evidence of transcription, which is particularly important for rare transcripts (e.g., and generally in species with few ESTs, such as rhesus and dog. For example, human DLX1-AS (antisense to DLX1) was Most samples exhibit higher correlation with cognate samples in other species than with samples in the same species (e.g., brain, liver, and testis). All clustering was performed hierarchically using Pearson correlation as a measure of distance and average linkage for grouping.
previously annotated only via its mouse orthologs (Supplemental Fig. 11B), and the let-7a-1 cluster in rhesus was previously annotated with a single EST (Fig. 3F). By extension, PolyA-seq results can be used to identify novel transcripts and can serve to complement RNA-seq data as well as data sets capturing 59 ends of mRNAs (Yamashita et al. 2011) and microRNAs (Chien et al. 2011). We observed a low level of polyadenylation in the body of some microRNA precursors (e.g., Supplemental Fig. 9B,C). Since transient RNA polyadenylation is reported to occur during RNA degradation (Slomovic et al. 2010), we hypothesized that PolyA-seq captures the 59 portion of primary transcripts following cleavage by DROSHA. We therefore examined the precise location of these polyA sites, and found that they lie immediately upstream of the mature transcripts (e.g., Supplemental Fig. 9C,D), consistent with transient polyadenylation of degradation precursors.
A major outstanding biological question is how alternative polyadenylation (APA) is regulated. In general, we observed similar proportions of usage between the dominant and alternative sites, even as absolute expression fluctuates (Supplemental Fig. 9). Even genes with splicing-dependent APA (e.g., Fig. 3C; Supplemental Fig. 13) tend to behave consistently across tissues. This supports the theory that APA is regulated by cis-acting regulatory effects such as sequence motifs and secondary structure within the transcripts and/or cis-acting DNA regulatory elements that impact transcription, which, in turn, affects 39 end processing (Di Giammartino et al. 2011). However, numerous studies have reported dynamic behavior of APA; 39UTR length can vary during development (Ji et al. 2009;Mangone et al. 2010), for example, and some 39UTRs become truncated in cancer samples . These and other observations support a second theory that it is the environment, namely, the concentrations of trans-acting factors, that regulate polyA site selection and usage. We believe that our five-species tissue atlas and PolyA-seq are uniquely positioned to help reveal the mechanisms underlying APA. To this end we envision the following future studies: 1. An array of approaches aimed at discovery of cis-acting regulatory elements, such as sequence motifs and nucleic acid structural elements, are possible. The atlas offers three dimensions of evidence for functional significance: quantitative levels of site usage, tissue specificity, and conservation of usage in other species. Motifs/structural elements that are enriched in many combinations of these three criteria (e.g., brain-specific motif present in humans and rhesus, but not in other mammals, that is associated with promoter-distal alternative polyA sites) will yield interesting biological insights into APA. 2. Correlation of APA with additional data, such as RNA-seq, ChIPseq, and DNase sensitivity maps to identify novel candidate trans-acting factors involved in regulating APA. Global genome binding maps of numerous transcription factors, histone modifiers, and RNA Pol II activities in many mammalian tissues are becoming available. A significant spatial correlation between binding affinity and a specific type of APA from the threedimensional sampling space of this atlas (see 1 above) may yield novel factors involved in APA. 3. Correlations of APA with gene expression to identify trans-acting APA regulators. Previous studies have revealed up-regulated levels of the core polyadenylation factors in iPS cells (Shi et al. 2009) and a general down-regulation in differentiated tissues, where 39UTRs tend to be longer (Ji et al. 2009), suggesting that this approach is worthwhile. An atlas of many tissues where both expression and polyA site usage are computed from the same data set offers a powerful resource to identify patterns of gene expression that correlate with APA (and possibly presence of cisacting factors). This may be a powerful approach to narrow the list of candidate factors thought to play a role in specific instances of APA. 4. Application of PolyA-seq to mouse knockouts or other genetic models, particularly systems where RNA-binding proteins are perturbed, could yield additional insight into APA. An example where this was done successfully was with NOVA2, which was ultimately shown to regulate APA (Licatalosi et al. 2008).
PolyA-seq as a method is also a practical approach for DGE and eliminates the need to account for transcript length when quantifying DGE, which is essential when using traditional RNAseq. PolyA-seq is thus better powered over RNA-seq to detect expression differences in short transcripts. Furthermore, since priming of reverse transcription is driven solely by adenines, biases in hybridization kinetics caused by transcript sequence differences are also averted. PolyA-seq can also accommodate paired-end reads, which could link polyA sites with upstream exons; we have done this successfully by simply changing one of the PCR primers to be compatible with paired-end flow cells (data not shown). One could also incorporate barcodes into the second-strand sequencing primer, which would further extend the scalability of the assay by enabling multiplexing. Application of PolyA-seq to disease states may reveal aberrant polyA-site usage, potentially leading to a novel class of biomarkers and candidates to explain disease etiology.

Tissues
FirstChoice Human Reference brain total RNA was obtained from Ambion (Applied Biosystems). UHR RNA was purchased from Stratagene Corp. Mouse and dog tissue total RNAs were purchased from BioChain. Total RNA from rat tissues and remaining human tissues were purchased from Zyagen. Rhesus tissue RNAs were provided by Merck & Co., Inc. PolyA+ RNA was isolated using the Dynabead mRNA Purification Kit according to the manufacturer's instructions (Invitrogen, cat # 610.06). Supplemental Table 2 contains a summary of all RNA samples used in this study.

PolyA-seq library preparation and sequencing
Two microliters of 0.1 uM-tailed dT primer T(10)VN was combined with 150 ng of polyA+ RNA in a final volume of 11 mL (see Supplemental File 1 for all primer sequences and detailed protocol). The primer-template mix was heated at 65°C for 5 min and chilled on ice before adding 9 mL of reverse transcription master mix (4 mL of 5x buffer, 2 mL of 10 mM dNTPs, 1 mL of 100 mM DTT, 1 mL of RNaseOUT, and 1 mL of SuperScript III enzyme). The 20-mL reversetranscription reaction was incubated at 40°C for 90 min, 70°C for 15 min, and cooled to 4°C. RNA template was degraded by adding 1 mL of RNase H (Invitrogen Corp.) and incubating at 37°C for 20 min, 75°C for 15 min, and cooling to 4°C. DNA was subsequently purified using the QIAquick PCR Purification kit and eluted with 65 mL of elution buffer (Qiagen, Inc.). For secondstrand synthesis, 60 mL of purified cDNA was added to 40 mL of Klenow master mix (12 mL of water, 10 mL of 10x NEBuffer 2, 5 mL of 10 mM dNTPs, 3 mL of 5 units/ul exo-Klenow fragment; M0212L, New England Biolabs, Inc.) and 10 mL of 10 uM-tailed random hexamer primer. The 100-mL reaction was incubated at 37°C for 30 min and cooled to 4°C. DNA was purified from the second-strand reaction by incubating with 1.8 volumes of Agencourt AMPure XP beads (Beckman Coulter) for 5 min, washing twice with 70% EtOH, and eluting with 50 mL of elution buffer. This step greatly reduced the number of clones with inserts <40 nt. For PCR amplification, 33 mL of purified second-strand synthesis reaction was combined with 17 mL of PCR master mix (10 mL of 5x Buffer 2, 1 mL of 25 mM MgCl2, 1 mL of 10 mM dNTPs, 2 mL of 10 uM forward primer, 2 mL of 10 uM reverse primer, 1 mL of ExpandPLUS enzyme; Roche Diagnostics Corp.). Samples were denatured for 2 min at 94°C, followed by two cycles of 94°C for 10 sec, 40°C for 2 min, 72°C for 1 min; 8 cycles of 94°C for 10 sec, 60°C for 30 sec, 72°C for 1 min; 15 cycles of 94°C for 15 sec, 60°C for 30 sec, 72°C for 1 min; and 72°C for 5 min to polish ends before cooling to 4°C. Double-stranded DNA was purified using AMPure XP beads as described above. Libraries were sequenced on an Illumina Genome Analyzer IIx using v4 reagents, and bases were called using the Illumina Analysis Pipeline (v1.4). PolyA-seq raw read sequences can be found at the NCBI Sequence Read Archive (submission SRA039286), while alignments and filtered sites are available at GEO (GSE30198). Note that PolyA-seq reads are antisense.

Read alignment
Reads were reverse-complemented, then aligned with SOAP2 (Li et al. 2009) to the respective genome and to the set of splice junctions generated from RefSeq genes (Wheeler et al. 2003), Ensembl genes, and UCSC KG. All transcript/gene databases used in this study were downloaded from UCSC on May 2, 2011. Predicted splice junctions from ESTs, GENSCAN, and N-SCAN predictions were also used for genomic regions where experimentally defined transcript evidence was not available. All possible splice junctions corresponding to one-and two-exon skipping events in gene/ transcript models above were represented; a minimum of 5 nt overlap per flanking exonic sequence was required for an alignment to be considered. Splice-junction alignments were converted into genomic coordinates, and reads that aligned uniquely were retained for further analysis. 3P-seq C. elegans reads (Mayr and Bartel 2009) were aligned with BLAT (Kent 2002), which is better able to cope with 39-A overhangs.

Filtering internal priming events and estimating FDR
To generate a map of genuine polyA sites we ran PolyA-seq on UHR RNA as above, but used a primer ending in T10 instead of T10VN, and sequenced to 76 bp in single-end mode. Reads were aligned using BLAT (Kent 2002). Supplemental Figure 1 outlines our procedure for constructing the filtering model. Cases with perfectly aligned 59 ends, but unaligned 39 adenines (three or more), were considered genuine polyA sites. Cases where the entire read aligned perfectly were considered potential internal priming events, excluding sites with an AAUAAA motif or its variants (Beaudoing et al. 2000) at À40 to À10 nt and/or sites overlapping a 39 end of a RefSeq transcript or EST (610 bases from the end). The two distributions of priming site base frequencies were used to generate a logodds model, where polyA Score ¼ log 10 pA and I are the genuine polyA site and internal-priming base frequencies, respectively. x i (A i ,C i ,G i ,T i ) denotes the base identity at position i within the downstream site being assessed, where A i , C i , G i , T i , are the frequencies of As, Cs, Gs, and Ts, at position i in the matrix. Thus, any 10-base sequence can be used to compute a polyA score, and performance on UHR T10VN data was excellent. The test set comprised 22,551 sites that agreed with genuine polyA sites above (positives) and 15,530 sites that agreed with internal priming events (negatives). At a score threshold of 3.0, we obtained a sensitivity of 85.6% and a specificity of 97.5% (1 -false-positive rate) on the test set.
Estimating the false-discovery rate (FDR; proportion of all predicted polyA sites that are incorrect) is more challenging than determining the false-positive rate (FPR; proportion of known negatives that are incorrectly called), since it requires knowledge about sites outside of the validated test set. Ideally, one would know where all the internal priming sites are and could measure FDR directly, since it would be exactly the same as FPR. Nonetheless, FDR can be estimated from the test set if the ratio of positives to negatives in the test set reflects the true underlying ratio in the complete set (i.e., the genome). In other words, if polyA sites and internal priming events are randomly sampled from the genome to construct the test set, FDR can be estimated as one minus the precision [i.e., the proportion of all predicted positives that are correct, or TP/(TP+FP)]. This lack of sampling bias is true for our test set; the ratio of polyA sites to internal priming sites was 1.45 (22,551/15,530) and the ratio of genuine polyA sites to internal priming stretches (measured from T10 reads) was 1.61 (0.92 million reads that are genuine polyA sites plus 5.7 million reads which are >98% As and likely derived from polyA tails, divided by 4.1 million reads that are likely internal priming events; see Supplemental Fig. 1); 1.61 is greater than 1.45, thus, we reduced the precision by this factor to estimate FDR (Supplemental Fig. 2D). At a polyA score threshold of 3.0, we estimate FDR to be 2.5%.

Comparison with known polyA sites
Coordinates of known transcript termini (Fig. 3B) were compiled from RefSeq, UCSC KG, and Ensembl transcript models (downloaded from UCSC) (Fujita et al. 2011), taking the 39-most site when transcript isoforms overlapped on the same genomic strand. Known 39UTR coordinates (Fig. 3D,E) were also compiled from RefSeq, UCSC KG, and Ensembl transcripts by collapsing overlapping isoforms (on the same strand) into a single 39UTR model that represented the union of genomic coordinates. A total of 122,215 UTRs from these transcript sources were collapsed into 27,175 UTR models, of which 20,873 were detected in at least one human sample (16,387 were detected on average in each tissue). Previously reported polyA sites (Fig. 3G) were compiled from PolyA-DB2 (Lee et al. 2007), dbEST (Boguski et al. 1993), GenBank mRNAs (Wheeler et al. 2003), UCSC Known Genes (Fujita et al. 2011), RefSeq genes (Wheeler et al. 2003), and Ensembl genes (Flicek et al. 2011). Coordinates were obtained from UCSC (Fujita et al. 2011) when available, or by aligning against the respective genomes using BLAT (Kent 2002) (using default settings and taking the best alignment).

Evolutionary conservation of sites
To minimize differences due to sequencing depth, we randomly removed mapped reads within each sample to match mouse testis (5.2 million mapped reads), so that all samples had the same number of input reads (except for mouse brain, which had 3.3 million in toto). We then translated site coordinates from mouse, rat, dog, and rhesus to human using the LiftOver program and corresponding chain files, all obtained from UCSC; strandedness was preserved during this step. Sites lying within 30 bp on the same strand were clustered, and read counts were summed for each sample within each cluster. Multispecies clusters were retained only if they contained a site detected in at least one sample in every species (n = 2590). For every pair of samples, we calculated the Pearson correlation coefficient of the natural log of read counts across all polyA site clusters. Samples were then clustered based on these correlations. We used LiftOver in a similar manner to determine the number of human polyA sites lying in regions conserved in all other species and the fraction supported by PolyA-seq data in the other species.