Conservation of alternative polyadenylation patterns in mammalian genes

Identifying conserved and non-conserved poly(A) sites

We performed a complete mapping of all 3' ESTs and full-length cDNAs onto the human and mouse genomes. After clustering EST and cDNA hits, potential poly(A) sites were identified based on several quality criteria including the presence of at least two ESTs/cDNA ending at site, reduced dangling ends in Blast matches, lack of potential internal priming tract in downstream genomic region and presence of a canonical or variant poly(A) signal near 3' end. We identified a total of 66,647 and 52,270 candidate poly(A) sites in the human and mouse genome, respectively which were then mapped to flanking Ensembl-annotated genes.

Alternative poly(A) site may be found in tandem in the same 3' exon, or in different 3' exons when associated to alternative splicing. We want to avoid the latter case, as poly(A) site usage may be dictated first by alternative splicing, which is itself conserved to some extent for specific genes in animal genomes [10, 11]. To avoid interference from alternative splicing, we only considered poly(A) sites located in the 3'-most exon, and following the 3'-most stop codon in case of alternative splice forms. This retained 27,654 poly(A) sites for 14,574 human genes and 25,987 poly(A) sites for 15,199 mouse genes. The resulting estimate of 1.8 poly(A) sites/gene is comparable to previous ones [4, 5] with a slight increase for mouse imputable to an expanded mouse EST database.

To compare poly(A) sites located on ortholog gene pairs in human and mouse, we aligned 3'UTR regions of all 14,481 ortholog pairs and we defined as "conserved" those human and mouse poly(A) sites displaying aligned poly(A) signals and EST/cDNA support in both species (Figure 1A). Poly(A) sites with EST/cDNA support whose poly(A) signals were not aligned were considered as non-conserved (Figure 1B), even when the cleavage site itself was properly aligned (Figure 1C). Multiple aggregated cleavage sites occurring after a single poly(A) signal (only ~1% of total poly(A) sites in our procedure) were discarded. We obtained a total of 4,807 conserved poly(A) sites and 33,458 non-conserved poly(A) sites. The conserved/non-conserved ratio is about 0.3 in either species. Among conserved poly(A) sites, 3711 were single sites and 1096 were multiple sites from 503 orthologous gene pairs. By our definition, 20% of human genes have a conserved poly(A) site and 2.5% of human genes have multiple conserved poly(A) sites. Figure 2 shows that conservation is higher in single poly(A) sites (33%) than in tandem poly(A) site (18%). This suggests that alternative poly(A) sites are evolutionally less conserved than single poly(A) sites. The complete list of human/mouse conserved poly(A) sites is presented in Additional file 1.

Number of sites and position in UTR

Tian et al. [4] have reported that poly(A) site configurations (single sites, tandem sites or sites located in different exons) tend to be conserved between human and mouse. The conservation of the number of tandem poly(A) sites in homologous gene, however, was not assessed. Figure 3 shows the distribution of poly(A) site numbers in human and mouse orthologous pairs, for genes having one or more sites, conserved or non-conserved. Numbers of poly(A) sites are significantly correlated in orthologs (P = 6.2 × 10^-260, χ² test). In other words, a human gene with multiple poly(A) sites tends to have multiple poly(A) sites in mouse too, suggesting that a selective mechanism acts on the number of alternative polyadenylation sites.

Do conserved tandem poly(A) sites show any positional preference in the 3'UTR region? We examined all tandem poly(A) sites, and classified them as "proximal" or "distal" according to their position relative to the stop codon. For genes with an odd number of sites, the site located in the central position was excluded (10% of sites overall). Figure 4 shows relative locations of conserved and non-conserved tandem poly(A) sites. In both human and mouse, conserved sites tend to occur more often in the proximal part of the UTR, while non-conserved sites tend to occur more often in the distal part.

Processing efficiency

We estimated the relative processing efficiency (RE) of poly(A) sites based on EST counts, normalized in such a way that the highest EST count of all poly(A) sites from the same gene had a value of one. This eliminates biases resulting from different EST coverage in different genes and between human and mouse (see Methods). Figure 5A compares the human/mouse correlation coefficients of RE in conserved poly(A) sites (r = 0.45, arrow) and in 10,000 sets of randomly selected non-conserved sites from orthologous genes (histogram). The relative efficiency of conserved poly(A) site is correlated between human and mouse, while that of non-conserved sites from orthologous genes is not.

Figure 5B shows the distribution of relative efficiencies in conserved and non-conserved poly(A) sites. If one considers as "major" any site with RE above 0.5, then conserved sites are more often of the major type (85%) than non-conserved sites (75%). This observation raises a question about the nature of conserved sites. Is the higher processing efficiency of conserved sites associated to some selective constraint, or to the presence of non-functional and/or misprediction among non-conserved sites? To answer this question, we subdivided non-conserved mouse sites according to their conservation or lack thereof in rat orthologous genes (Figure 6). Poly(A) sites conserved between mouse and rat behave in the same fashion as poly(A) sites conserved between mouse and human, i.e. with a predominance of "major" sites. This suggests that the higher efficiency of conserved human/mouse sites is due firstly to their status of biologically functional sites, rather than to a property of ancient conserved sites.

Tissue specificity

We then analyzed the tissue specificity of poly(A) sites based on the eVOC expression ontology system [12]. Version 2.6 of eVOC maps each of the 9,478 human EST libraries to a formalized tissue description. As a mouse version of eVOC was not available at the time of the study, we further mapped 556 mouse EST libraries using the same formal description system (See Materials and Methods). We obtained for each poly(A) site the number of different tissues in which the site is observed, among the 12 possible top-level eVOC tissue categories. Since tissue counts are highly dependent on EST coverage, we normalized tissue counts versus an expected number of tissues obtained from a random EST sample of same size. Our measure of tissue specificity is the log ratio of observed vs. expected number of tissues. About 10% of tandem poly(A) sites have a tissue specificity below -0.5 (high specificity) and <3% above 0.5 (low specificity). Non-conserved poly(A) sites showed no correlation in tissue specificity (Fig 7A, histogram), while conserved sites had weakly correlated tissue specificities (r = 0.10, Fig 7A, arrow). This observed difference is not significant based on a T-test performed after Fisher's z-transformation of the r value (P < 0.5).

To circumvent possible gene-level expression biases, we measured a "relative tissue specificity" (RTS), by assigning a value of 1 to the poly(A) site with the broadest tissue distribution in a gene. Each gene thus has at least one site with RTS = 1. The distribution of other sites is shown in Figure 7B. Interestingly, very few sites have a RTS below 0.5. The fact that most sites have a RTS close to 1, that is close to the broadest possible tissue distribution for this gene, means variations in tissue specificity between successive poly(A) sites in a gene are generally limited. We used the median RTS to distinguish "broad" from "narrow" sites (we preferred these terms over "constitutive" and "specific" since these would suggest an absolute usage level, while here we only measure relative usage). Sites with an RTS above median (0.90 for human, 0.88 for mouse) are said to display a "broad" tissue distribution while other sites are said to display a "narrow" tissue distribution. Based on this definition, broad and narrow tissue distributions are equally frequent among conserved and non-conserved sites (Figure 7B). Broad and narrow tissue distributions are also equally distributed among major and minor sites (Additional file 2).

Spatial preferences vs. efficiency or specificity

We examined the relationship between the spatial organization of tandem poly(A) sites and site efficiency or specificity. Poly(A) sites were classified as major/minor or narrow/broad as described above, and their spatial organization was observed in genes containing two or more tandem conserved sites. A "conserved usage pattern" was recorded when successive tandem sites had the same efficiency and/or specificity pattern in human and mouse orthologs (Figure 8A). We observed that 53% of genes with tandem conserved poly(A) sites had a conserved efficiency pattern (227 gene pairs), significantly higher than expected by chance (146 gene pairs, binomial distribution P = 9.7 × 10^-28). Comparing expected and observed values, 81 genes would have a poly(A) site efficiency pattern under selection. This tendency is not observed for tissue specificity: 133 gene pairs have a conserved tissue specificity pattern vs. 141 expected by chance.

In Figure 8B, we focus on genes containing at least one conserved major site, as a surrogate for sites existing prior to human-mouse divergence. We then observe how such sites are associated to flanking conserved (black) or non-conserved (gray) sites, using flanking non-conserved sites as surrogates for emerging sites. Interestingly, while emerging minor sites are as frequent on the 5' or 3' side of existing major sites (top row), pairs of conserved sites (bottom row) are twice more often 5'minor/3'-major than 5'-major/3'-minor. This suggests that selection of alternative poly(A) sites favors the pattern 5'-minor/3'-major over the pattern 5'-major/3'-minor.

Differentially processed Poly(A) sites

We define here as differentially processed those alternative poly(A) sites with a significantly biased usage in any tissue class, as compared with other poly(A) sites from the same gene. Differential usage was measured using a Fisher test as previously reported [6]. Using a Bonferroni correction for multiple testing, five conserved and 84 non-conserved sites are differentially processed in human (54 and 369, respectively, in mouse).

We examined the relationship between differential processing and site efficiency or tissue specificity. Consistently with a recent study of tissue-specific polyadenylation [8], minor sites are more often differentially processed than major sites (Figure 9A). Although this study did not identify tissue biases in major sites, we do observe a few occurrences (9 in human, 64 in mouse) of major sites with differential processing. Differential processing is also much more frequent in non-conserved than in conserved sites (Figure 9B).

Expectedly, there is a high correlation between differential usage and the "narrow" or "broad" status of a site. Differentially processed sites are about three times more often of the narrow type than of the broad type (data not shown). Although counterintuitive, some poly(A) sites can be at the same time differentially processed and of broader usage, because our specificity measure is relative and always classifies as broad the site with the broadest tissue usage, even when usage is restricted to a single tissue.

A list of differentially processed, conserved poly(A) sites is presented in Additional file 3. Differentially processed sites are observed in all tissue classes (Figure 10). The apparent overrepresentation of urogenital and nervous systems is not significant when EST library coverage is taken into account. EST coverage is not sufficient either to provide interspecies confirmation of tissue biases. No conserved site is found differentially processed in both human and mouse after Bonferroni correction.

Conserved sequence motifs around conserved sites

As our criteria for poly(A) site conservation imply a correct alignment of poly(A) signals, we suspected that conserved sequences around poly(A) signals could also contribute to poly(A) site conservation. This region is known to contain elements such as the USE (upstream sequence element) and DSE (downstream sequence element), two U-rich elements involved in the control of poly(A) site efficiency [13–15], as well as a number of potential regulatory motifs of unknown function [7]. A possible explanation for proper signal alignment and increased cleavage efficiency at conserved poly(A) sites could be related to the occurrence of such control elements in both human and mouse orthologs. Although downstream regions appear slightly more U-rich in conserved sites than in non-conserved sites (see Additional file 4), indicative of stronger DSE elements [15], we could not find overrepresented sequence motifs occurring in more than a few conserved sites. Therefore, there is no widespread cis-regulatory element that would explain poly(A) site conservation.

Poly(A) site prediction

EST sequences were obtained from dbEST v. 01/06/05 (6,009,051 human, 4,314,509 mouse and 623,741 rat ESTs). Full-length cDNA sequences were obtained from H-Inv 1.8 (41,118 sequences) [18] and FANTOM 2.01 (60,770 sequences) [19]. ESTs annotated as 3' were extracted (2,029,534 human, 1,864,874 mouse and 293,597 rat ESTs) and trailing poly(A) or poly(T) sequences of 5 nt or more were removed, with one mismatch (non-A or non-T) allowed for polyA/T tails of 10 or more. Both 3' EST and cDNA sequences were aligned to the repeat-masked human genome v27.35a.1, mouse genome v27.33c.1 and rat genome v27.3e.1 using the Megablast program [20]. We did not use a specific exon junction mapping software, since we were only interested in the terminal part of the 3' exon. All hits presenting at least 95% identity with the genomic sequence were retained (hit size >28 nt at default E-value). Partial hits flanking a repeat masked region of the genome were then realigned to the locally unmasked region. Hits with 95% identity after this step were retained. Clusters were formed with ESTs having either their 5' or 3' extremities falling within a 10 nt distance. ESTs were not oriented at this stage. Each cluster was analyzed using a sliding window to locate the most likely cleavage site, defined as the position where the window contains the most EST/cDNA ends. The following filters where then applied:

(i) Dangling ends: discard hits with more than 5 unmatched nt at cleavage site

(ii) Internal priming: discard cleavage sites flanked by A-rich region (at least 9 As out of 10 nt) in the 50 nt downstream genomic sequence

(iii) Poly(A) signal: retain only cleavage sites where 30 nt upstream genomic sequence contains one of the 11 variant poly(A) signal identified in our previous study [2]: AATAAA, ATTAAA, AGTAAA, TATAAA, CATAAA, GATAAA, AATATA, AATACA, AATAGA, AATGAA, ACTAAA

Only those cleavage sites passing the filters and supported by at least two ESTs/cDNAs were retained as predicted poly(A) sites. From the starting EST/cDNA datasets we finally retained 718,927 ESTs and 19,626 full length cDNAs to identify 66,647 different poly(A) sites in human, 739,259 ESTs and 23,069 full length cDNAs to identify 52,232 different poly(A) sites in mouse and 119,946 ESTs to identify 27,494 different poly(A) sites in rat.

Assignment of tandem Poly(A) signal sites

Poly(A) sites were assigned to transcript sequences taken from Ensembl 27.35a.1 [21]. If the poly(A) site lied within one or more annotated transcripts, downstream of the end of translation, the site was affected to each of these transcripts. If the poly(A) site lied upstream of the end of translation, then it was considered as "in CDS" and is not used for analysis. If the poly(A) site did not map to any annotated transcript, it was affected to the nearest 5' transcript. Poly(A) signals were assigned to their respective poly(A) sites by taking the signal that was closest to the 5'-most poly(A) site in each cluster. Only poly(A) sites mapping to the 3'-most exon of an Ensembl gene or its genomic downstream region up to 10 kb were considered further.

Assignment of conserved Poly(A) sites

Ortholog human/mouse (or mouse/rat) gene pairs were obtained from EnsMart [22]. All genes with paralogs were omitted from the analysis. 3'UTR regions assigned in Ensembl including up to 10 kb downstream genomic sequence of all transcripts were aligned by ClustalW with default parameter. Predicted Poly(A) sites were then defined as conserved if they were within a distance of 30 bp of a properly aligned poly(A) signal and had EST-support in both human and mouse. In the case where multiple poly(A) signals were associated to a single cleavage site, the signal closest to the cleavage site was used for the analysis. Multiple cleavage sites that were associated to the same poly(A) signal were omitted.

Mouse eVOC ontology mapping

Anatomical terms in mouse cDNA libraries were mapped to anatomical systems from the eVOC ontology ver 2.6 [12]. Terms that matched exactly to human eVOC terms were kept and all mixed terms (e.g. lung and colon) were manually checked to see if each component could match to a human eVOC term. Any terms that could not be directly mapped were classified as "unclassifiable". Finally all libraries were assigned to 12 anatomical systems: anatomical site, cardiovascular, respiratory, hematological, lymphoreticular, alimentary, urogenital, endocrine, musculoskeletal, dermal, nervous and unclassifiable.

Correlation analysis

χ² test, t-test, calculation of correlation coefficient and Fisher's z-transformation were performed using Microsoft Excel 2002. Distributions of observed and random distributed values were calculated using dedicated Perl scripts.

Efficiency and tissue specificity

All EST counts were performed after discarding EST libraries annotated as normalized in the dbEST database (3% of overall human and mouse libraries). The relative efficiency R of a poly(A) site was calculated as a ratio of number of ESTs,

$R_i=\frac{n_{X,i}}{n_{X,\max }}$

where n_x,iis the number of ESTs of poly(A) site i within gene X, and n_X,maxis maximum number of ESTs of any tandem poly(A) site within gene X.

Sites with a ratio higher than 0.5 were defined as "major" (high efficiency), while other sites were defined as "minor" (low efficiency).

The specificity of a poly(A) site was defined as the number of different expression systems in which this site was utilized, normalized by the number of supporting ESTs as follows. For a poly(A) site supported by N ESTs, we calculated an expected number of expression systems as the average number of expression systems obtained in 10,000 random sets of N ESTs sampled from the complete set of ESTs mapping any poly(A) site. Specificity S was then calculated as the log-ratio of the number of expression systems:

$S=\log \left(\frac{T}{T_{sim[n]}}\right)$

where T is observed number of different expression system in each poly(A) site, T _{sim [n]} is the simulated number of different expression systems for n ESTs.

For genes with tandem poly(A) sites, a relative tissue specificity U was calculated for each site, as the ratio of tissue specificities:

$U_i=\frac{S_{X,i}-S_{\min }}{S_{X,\max }-S_{\min }}$

where S_X,iis tissue specificity of each tandem poly(A) site i within gene X, S_X,maxis maximum tissue specificity of tandem poly(A) site within gene X, and S _min is minimum value of S in each species. We adjusted the interval of relative tissue specificity to 0–1; 0 being the most specific site and 1 being the least specific. Median ratio was 0.90 for human, 0.88 for mouse. Sites with a ratio higher than these values were defined as "broad" while others were defined as "narrow".

A usage pattern was defined as the sequence of relative efficiencies or relative tissue specificities for all tandem poly(A) sites in a gene. For each usage pattern, the expected value E of the number of gene pairs that could randomly bare this pattern was calculated from a random combination of human and mouse gene pairs with conserved tandem poly(A) sites as described below. For a number of conserved poly(A) sites j, the list of relative usage patterns L [j, k] was defined by

L [j, k] = {[r₁, r₂,...,r_j]₁, [r₁, r₂,...,r_j]₂, ..., [r₁, r₂,..., r_j]_k}

where r_j = {major or minor} or {broad or narrow}. For a number N of genes and the number of conserved poly(A) sites j, a probability P for human and mouse poly(A) sites having pattern L [j, k] was defined as:

$P_{\text{human}}=\frac{N_{L[j,k],human}}{N_{total}},P_{\text{mouse}}=\frac{N_{L[j,k],mouse}}{N_{total}}$

Then the expected value E for the maximal value of k (k_max = 2^j-1) is :

$E={\displaystyle \sum _{j=1}^{j_{\max }}\left\{{\displaystyle \sum _{k=1}^{k_{\max }}\left(N_j\ast P_{human}\ast P_{mouse}\right)}\right\}}$

Expected frequencies of genes with narrow usage patterns, based on a binomial distribution, were calculated by Microsoft Excel 2002.

Differential use of poly(A) sites

To identify differentially processed poly(A) sites, Fisher's tests were performed on the distribution of the number of supporting ESTs from each expression system against all other systems for each poly(A) site as previously described [6]. A Bonferroni correction for multiple testing was applied. Poly(A) sites supported only by ESTs from pooled tissue libraries were omitted.