Two to four percent of RefSeq genes are nonsense-mediated decay candidates and not explained as rare alleles or sequencing artefacts

By reference to the Mammalian Gene Collection full length cDNAs [49, 50] and dbSNP [51] data (Table S1 in Additional file 1), we find that the NMD candidates have the same quality support as other genes, indicating that these candidates are not the result of sequence artefacts or rare alleles. We also found that human NMD candidates here are enriched for NMD potential targets, previously determined in the study by Mendell et al [10] (Table S2 in Additional file 1). As shown in Table 1, NMD candidates generally encode shorter proteins and longer introns and UTRs compared with non-NMD genes.

Few genes are nonsense-mediated decay candidates in both mouse and human

As it is likely that some NMD transcripts were missed from the analysis due to strict data selection, we repeated the analysis using all the RefSeq mRNAs, including predicted mRNAs (with prefix 'XM_'). The number of NMD candidates in both genomes increased as expected in this second round. However, the intersection between NMD orthologs was still small (see Table S3 in Additional file 1). The result is consistent with a previous comparison of human, fruitfly D. melanogaster and S. cerevisiae NMD candidates [4, 9–12]. We can also show that low rates of NMD conservation are not consistent with the normal rates of stop-codon turnover (see Table S4 in Additional file 1). Our finding shows that, even when the mechanism of PTC regulation is not a variable, deterministic regulation by NMD is generally not selectively favoured over the long term or is not the correct explanation for most PTCs.

Nonsense-mediated decay exons tend not to be multiples of three long

If an exon is included by noisy splicing, pressure not to be translated should be higher if the exon is not a multiple of three long than if it is a multiple of three long. The inclusion of such exons will induce a frame shift if translated, which will change the encoded amino acids downstream of included exons. If so, this will result in proteins that are at best costly to make and non-functional, and at worst are toxic to cells. In contrast, inclusion of an exon that is a multiple of three will typically result in a small peptide insert but need not disrupt the overall function of that protein. The costs of noisy splicing can be reduced to some extent by degradation of these noisy transcripts employing mRNA decay systems, such as by NMD. As exons that are not a multiple of three impose a higher cost, we expect stronger selection for PTCs in exons that are not multiples of three compared with those that do not induce frameshifts, assuming an equal rate of mis-splicing. The regulation model makes no prediction of bias.

Nonsense-mediated decay-specific exons tend to be in the low inclusion category and newly created

According to the noisy splice model, the NMD transcript is an alternatively spliced unwanted transcript. The regulation model does not require the NMD transcripts to be alternatively spliced transcripts, nor, if they are, need they be the minority form (the transcript isoform that constitutes a small fraction (less than one third) of transcripts from the same gene). Are then NMD genes more likely to be alternatively spliced than random genes and are the NMD transcript isoforms rare? To investigate this we mapped our gene lists to Ensembl genes with BioMart [53], and extracted the splicing isoform information from the ASD database [54, 55].

As expected from the noisy splicing model we found that a larger proportion of NMD-specific exons (48%) were spliced in minor form (included in less than one third of transcripts transcribed from this gene) compared with non-NMD-single exons (26.1%) (Figure 2) (One-sided Kolmogorov-Smirnov test, P = 0.02443). This is consistent with previous reports that most putative NMD transcripts are expressed in low abundance across examined tissues [41] and with the finding that both in humans and mice the minor transcript has a PTC more commonly than the major form (11.1% versus 3.7%) [44].

Are NMD exons ancient or new? To explore this question, we examined the exon creation and loss events for each NMD-specific cassette exon. To obtain this information, we started by mapping our exons to those from the ASAP2 and the VEEDB databases [56, 57]. The VEEDB database provides exon conservation information for each given exon based on splice site conservation, this being extracted from a 17-vertebrate UCSC multi-genome alignment [57]. Using this exon conservation data, we could determine whether a given human exon is conserved or absent in mouse and dog (outgroup). Unfortunately, only a small portion of exons are mapped to the VEEDB database [57]. As shown in Table 4, 9 out of 35 NMD-specific exons were created after the human-mouse split. This proportion (25.7%) is significantly higher than that of non-NMD-single class exons (11.9%) (Fisher's exact test, P = 0.0315) (see Additional file 2 for human NMD-specific exons conserved in mice). When we compared humans against Rhesus macaques (Macaca mulatta) with mice as outgroup, the difference is more significant (Additional file 1, Table S5, NMD-specific 21.2%, non-NMD-single 4.5%, P = 0.001139). These findings are consistent with the fact that NMD-inducing exons are often not conserved among species [41]. The exon loss rates are much smaller and there was no difference between NMD and non-NMD exons.

An association between new and alternatively spliced exons can probably also account for the rapid turnover of genes subject to NMD. Two types of species-specific alternative splicing events can be defined [58, 59]. One type, referred to as 'species-specific alternative splicing of conserved exons', is represented by a conserved exon that is alternatively spliced in one species but constitutively spliced in the other species. The other type, referred to as 'genome-specific alternative splicing', is represented by an alternative exon in one species which is not detectable in the ortholog of the other species (see Figure Three in [41] for diagrams). More than 41% of species-specific alternative splicing events of conserved exons and 61% genome-specific alternative splicing events had the potential to trigger NMD, while the NMD-inducing rate of conserved alternative splicing events between human and mice was much lower (< 31%) [41]. No matter which form it is, both types cause NMD to occur more often in only one of the two species than in both species, and can hence explain to a considerable extent NMD status divergence.

Where an orthologous exon can be found, the nonsense-mediated decay-specific exon is or was under purifying selection

Comparing alternative cassette exons with other parts of the same gene may be an error-prone test as cassette exons may generally have weaker evolutionary pressures due to their exclusions in some splicing isoforms. Indeed, we observe a higher K _a/K _s ratio in non-NMD-single class of exons than that in other regions (Table 6). To exclude the effect of this, we compared the NMD-specific K _a/K _s ratio against non-NMD-single exons. A higher K _a/K _s ratio for the NMD-specific class was still observed (Table 6, 0.3423 versus 0.2554).

While these results are consistent with the regulated expression model, fuller interpretation of the results is non-trivial. First, the ability to detect orthologous exons predisposes to finding exons functioning in regulation. Second, even if the spurious transcript model is correct for the lineage with the stop, the null for K _a/K _s is not 1. Crucial is when the exon was subject to NMD and what form of selection operated prior to this. If the exon was subject to NMD very recently, then most of the evolutionary history of the exon down the human lineage was not evolving in response to the presence of the stop. Only if the exon was always spurious, unlikely given that this exon is found in multiple distant taxa, would K _a/K _s = 1 be expected for the lineage in question. In short, by virtue of the fact that we can find distant orthologous exons, we are almost certainly biasing the data set to exons that have been and possibly still are functional. Given that K _a/K _s < 1 we can be confident that for some of the time the exon has not been spurious.

Nonsense-mediated decay candidates are fast evolving but not hot-spots for adaptive evolution

Given that NMD candidates evolve faster in their own lineage than orthologs (relative rate test, Table S7 in Additional file 1), possibly owing to reduced selective constraints, it is tempting to suppose that NMD genes and exons are potentially given much more freedom to roam sequence space than exons of, for example, house-keeping genes. Might this predispose NMD genes to be hot-spots for adaptive evolution?

Under-representation of some functional classes of genes in the nonsense-mediated decay set is consistent with noisy splicing

To examine the skew we employed tools of PANTHER database [65, 66]. The biological processes with Bonferroni-corrected P values < 0.1 in either species are listed in Table 9. About half of the NMD candidates were classified as Biological process unclassified. This was the only set showing over-representation with the NMD class. The only other classes showing significant deviation from expected showed under-representation, these being Developmental processes, Cell-surface receptor-mediated signal transduction, Sensory perception, Signal transduction, Chemosensory perception, and Olfaction. The functional distributions of NMD candidates in humans and mice were quite similar. Of the top six most significant GO terms from either species, five (Biological process unclassified, Developmental processes, Sensory perception, Signal transduction, Chemosensory perception) also appeared in the list of the most significant in the other species (Table 9).

To further determine whether the divergence of orthology in NMD candidates also led to functional divergence, we compared the functional distributions of human and mouse NMD candidates using the FatiGO web tool [67], which is able to detect particular GO terms for which the two lists of genes have different proportions of genes. For feasibility, the GO terms for mouse genes were deduced from the corresponding human orthologs. No significant GO terms were detected at any level (GO levels 3 to 9; see Additional file 4). To exclude the effect of orthologs that are NMD candidates in both species, we repeated the analysis after removing these orthologs from either or both species. As before there was no GO term showing a significant difference between mice and humans in the regularity of NMD (data not shown). This indicates that there is no functional class in which there are significantly more or fewer human NMD candidates compared against mouse NMD candidates.

These results suggest that NMD targets different genes in the two species, but ones largely in the same functional categories (Table 9 and Additional file 4). While superficially this looks like evidence for regulated splicing, if alternatively spliced genes are more prone to incorrect splicing, we expect, under the noisy splicing model, that the NMD-under-represented functional classes will have fewer alternatively spliced genes than other classes.

Data collection

We downloaded the sequences and annotations for human and mice from the NCBI RefSeq [46] database (Build 36.1) in January 2007. To improve the confidence of NMD candidate identification, we only extracted the transcripts with initial letters 'NM_'. Based on these annotations, for each transcript, we calculated the distance from the stop codon to the exon-exon junction closest to the 3' end. According to the NMD rule [3, 4, 16], we classified transcripts as NMD candidates if the distance was > 55 nucleotide bps. Then, we defined a gene as an NMD candidate if it had at least one NMD candidate transcript. All the remaining genes in the genome were classified as non-NMD candidates.

The UTR length, average intron length and protein length for each transcript were also calculated or extracted from the annotations. For each gene with multiple transcript variants, we collapsed these parameters into one by choosing the splicing form with the longest protein and calculating the means of transcript isoforms.

Exon type classification

We started from 3,362 genes with at least two RefSeq mRNAs in our dataset, which included 252 NMD and 3,110 Non-NMD candidates. We only considered the coding exons in each transcript and excluded the two marginal 5' and 3' exons within each gene due to their general incompleteness. For NMD candidates, we searched the exon isoforms only observed in NMD transcripts and defined these exons as NMD-specific exons, and defined the remaining exons in NMD candidates as NMD-non-specific exons. For non-NMD candidates, we classified exons as non-NMD-single and non-NMD-multiple. The former were observed in only one splicing transcript isoform for a given gene and the latter were observed in at least two different splicing transcript isoforms.

Mapping exons to ASAP2 database

Since the ASAP2 database gives the positions of exons on human chromosomes [71] of NCBI build 35.1 [46, 72, 73], equivalent to UCSC hg17 [74–76], we mapped the human exon set of RefSeq mRNAs to the ASAP2 database [56] as follows: first, we converted the exon positions on reference sequence contigs into those on chromosomes (NCBI build 36.1, UCSC hg18) using a Perl script. Second, we used the UCSC [76] liftOver tool to convert these positions into those on human NCBI build 35.1 (UCSC hg17). Finally, we compared these positions with those in the ASAP2 database and retrieved the ASAP2 exons that exactly matched the RefSeq exon set. 95,029 of 209,222 RefSeq exons can be uniquely mapped to the ASAP2 database. Based on these matched exons, we can easily obtain the splicing state, inclusion levels, exon creation/loss from the ASAP2 and the VEEDB database tables [56, 57].

Construction of exon alignments for nonsense-mediated decay -specific exons and calculation of K _a/K _s

Given the NMD-specific cassette exon lists above, we extracted the corresponding regions from human-mouse-dog CDS alignments (built by a Clustal W, version 1.83 [77], see Additional file 1 for details) and concatenated them together. The remained regions in these alignments were also concatenated. These concatenations were inputted separately into PAML [60] package to calculate K _a/K _s ratio for each lineage under the free ratio model [61]. The alignments for non-NMD-single exons and remaining parts were similarly extracted and inputted into PAML for K_a/K_s calculations.

To see if the NMD-specific exons were under neutral evolution, we fixed the K _a/K _s at one in the human NMD lineage under free ratio model [61], and compared this with a more general model (with K _a/K _s free) to test if the neutral model could be rejected using likelihood ratio tests (calculated in R [78]).

Gene ontology analysis

Comparisons of functional distributions of NMD candidates between human and mouse were carried out using the FatiGO program [67]. FatiGO implements the nested inclusive analysis, in which the test is done recursively until the deepest level in which significance is obtained and only this last level is reported. In this way both variables, the efficiency of the test and the highest precision in the term found, are optimized. The program computes a Fisher's two-tail exact test to statistically define over- or under-represented terms between two lists of genes, and the original P values are corrected by a false discovery rate approach [79].

The detection of over- or under- represented functional entries for NMD candidates was done based on the PANTHER database [65, 66]. The NMD candidate list was compared with the all the genes used and the P values were determined using a binomial test for each functional category. The original P values were adjusted using a modified Bonferroni correction method, which accounted for the nesting relationship among GO terms at different levels to avoid too conservative corrections.