(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



Transcript diversity reflects deleterious RNA processing errors shaped by population size in metazoans [1]

['Kai Mi', 'Bio-X Institutes', 'Key Laboratory For The Genetics Of Developmental', 'Neuropsychiatric Disorders', 'Ministry Of Education', 'Shanghai Jiao Tong University', 'Shanghai', 'Lili Guan', 'Department Of Otorhinolaryngology Head', 'Neck Surgery']

Date: 2026-03 

In eukaryotes, alternative transcription initiation (ATI), alternative splicing (AS), and alternative polyadenylation (APA) result in multiple different transcripts per gene, but the biological significance of the transcript diversity produced remains controversial. Some suggested that this diversity is adaptive, while others contended that it is largely deleterious and arises from molecular errors in transcription and RNA processing. The error hypothesis makes a distinct prediction that is not expected under the adaptive hypothesis: transcript diversity declines with the effective population size (N e ) of the species because natural selection minimizing errors is more effective under larger N e . By analyzing 166 transcriptomes from 75 metazoans, we report that transcript diversity measured by the percentage uses of minor ATI, AS, and APA sites decreases with N e or its proxies. This observation supports the error hypothesis and suggests that metazoan transcript diversity is largely deleterious.

Funding: This study was supported by National Science and Technology Innovation 2030 Major Projects for ‘Brain Science and Brain-Inspired Research’ (2022ZD0214400 to CX), National Natural Science Foundation of China (32270704 and 32472630 to CX), Medical-Engineering Crossover Fund of Shanghai Jiao Tong University (YG2025QNB51 to CX), and the U.S. National Institutes of Health (R35GM139484 to JZ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

In the present work, by quantifying ATI, APA, and AS in 166 transcriptomes, we compared levels of transcript diversity among 75 metazoan species spanning a wide range of N e . We found that transcript diversity generally decreases with N e or its proxies, strengthening the previous support of the error hypothesis for AS [ 36 ] and providing new evidence for this hypothesis for ATI and APA.

Nonetheless, the level of transcript diversity has not been extensively compared across species. Such a comparison is useful for differentiating between the adaptive and error hypotheses, because the two hypotheses make distinct predictions about the relationship between the level of transcript diversity in a species and the effective population size (N e ) of the species. Specifically, under the error hypothesis, transcript diversity is due to deleterious molecular error, so is disfavored and lowered by natural selection. Because the efficacy of natural selection increases with N e , we expect transcript diversity to decline with N e [ 41 , 42 ]. Under the adaptive hypothesis, however, transcript diversity is beneficial so may be selectively elevated. Under this scenario, transcript diversity is expected to increase with N e . However, one could also argue that, under the adaptive hypothesis, the optimal level of transcript diversity in a species depends on the specific condition and environment of the species; as a result, no prediction can be made regarding the relationship between the level of transcript diversity and N e . At any rate, a negative correlation between N e and transcript diversity is predicted by the error hypothesis but is not expected under the adaptive hypothesis. Indeed, this prediction was validated by a comparative analysis of AS across 53 species [ 36 ].

Despite the universality of ATI, APA, and AS in eukaryotes, the biological significance of the created transcript diversity is debated. Some case studies suggested that different RNA isoforms are functionally distinct. For instance, the human Lef1 gene [ 27 ], mouse Ighm gene [ 28 ], and Drosophila Sxl, Tra, and Dsx genes [ 29 ] have functionally distinct RNA isoforms (and corresponding protein isoforms) produced by ATI, APA, and AS, respectively. Such examples led to the hypothesis that transcript/proteome diversity is generally adaptive and that ATI, APA, and AS are widely used, regulated mechanisms to expand transcript/proteome diversity [ 3 , 4 , 30 , 31 ]. However, a competing hypothesis known as the error hypothesis has also been suggested [ 32 – 38 ]. The error hypothesis contends that transcription and RNA processing are error-prone; consequently, the vast majority of the observed transcript diversity reflects molecular errors that not only lower the number of functional molecules and waste energy but may also create cytotoxicity [ 39 ]. Several lines of evidence support the error hypothesis [ 32 – 36 , 39 , 40 ]. For example, the error hypothesis predicts that transcript diversity is lower in relatively highly expressed genes than in relatively lowly expressed genes because of stronger selection minimizing error rates acting on highly than lowly expressed genes [ 39 ]. Empirical data indeed support this prediction [ 39 ].

In eukaryotes, alternative transcription initiation (ATI) [ 1 ] and alternative polyadenylation (APA) [ 2 ] can respectively vary the beginning and end of a transcript produced from a gene, while alternative splicing (AS) can generate different RNA isoforms by selective inclusion or exclusion of exons in mRNA processing [ 3 ]. As a result, multiple different RNA transcripts (isoforms) are often produced from a single eukaryotic gene [ 4 , 5 ], generating transcript diversity. These transcripts may vary in their coding sequence, untranslated regions, and/or other regulatory elements [ 1 – 3 ]. ATI, APA, and AS are common phenomena in various eukaryotes such as fungi [ 6 – 8 ], plants [ 9 – 11 ], and animals [ 12 – 14 ]. For example, in humans, > 70% of genes exhibit APA [ 13 ], > 50% of genes display ATI [ 15 ], and >95% of multi-exon genes show AS [ 16 ], resulting in >170,000 transcripts recorded for ~20,000 human protein-coding genes (ENSEMBL genome reference consortium human build 38; GRCh38). ATI, APA, and AS may vary among tissues [ 17 – 19 ], across developmental stages [ 19 – 21 ], and during cell differentiation [ 22 – 24 ], and can contribute to disease [ 1 , 25 , 26 ].

(A–D) Relationship between transcript diversity caused by AS and N e (A) , life span (B) , body length (C) , or ω (D) in the brain. P-values from Spearman’s correlation and PGLS are shown. N represents the number of species included. (E) Correlation between transcript diversity caused by AS and N e (or proxies) in three tissues. BUSCO genes are used in all panels. The data underlying this Figure can be found in https://doi.org/10.5281/zenodo.18514977 .

Next, we calculated the mean transcript diversity (caused by AS) per gene for each species among BUSCO genes using data from the brain tissue. We found that this quantity significantly decreases with N e ( Fig 4A ), but increases with life span ( Fig 4B ), body length ( Fig 4C ), and ω ( Fig 4D ). Similar results were observed for the ovary and testis ( Fig 4E ). These patterns remain qualitatively unchanged ( S3B Fig ) when all protein-coding genes were analyzed. Thus, consistent with a previous study [ 36 ], our across-species comparison of AS supports the error hypothesis.

The error hypothesis predicts that the transcript diversity of a gene caused by AS should decrease with the expression level of the gene, because more highly expressed genes are subject to stronger selection against splicing error [ 64 ]. Indeed, we observed negative correlations in 158 of 166 samples ( S3A Fig ). These results confirm the previous findings from a limited number of species [ 34 ] and suggest that the error hypothesis of AS is broadly supported in animals.

Previous studies have shown that AS is noisy [ 62 ] and mostly nonadaptive [ 34 ], consistent with the error hypothesis. To compare transcript diversity caused by AS across species, we assembled the transcripts for a gene and quantified the expression level of each transcript of the gene using StringTie, a widely-used tool outperforming others in the accuracy of both assembly and expression level measurement [ 63 ]. We then computed for each gene its transcript diversity due to AS by dividing the total splicing amount of all minor RNA splicing isoforms by the total splicing amount of all RNA splicing isoforms of the gene. Here, the splicing amount of an RNA splicing isoform is the total number of reads covering all splicing junctions of the isoform.

(A–D) Relationship between transcript diversity caused by APA and N e (A) , life span (B) , body length (C) , or ω (D) in the brain. P-values from Spearman’s correlation and PGLS are shown. N represents the number of species included. (E) Correlation between transcript diversity caused by ATI and N e (or proxies) in three tissues. BUSCO genes are used in all panels. The data underlying this Figure can be found in https://doi.org/10.5281/zenodo.18514977 .

Next, we calculated the average total percentage usage of minor ATI sites per gene for the BUSCO genes of a species using RNA-seq data and correlated it with N e across species. Indeed, a significant, negative correlation was observed across 19 species (ρ = −0.74, P = 2.7 × 10 −4 ; P PGLS = 6.8 × 10 −3 ; Fig 3A ). We similarly observed a positive correlation when N e is replaced with life span (ρ = 0.66, P = 1.3 × 10 −5 ; P PGLS = 3.1 × 10 −2 ; Fig 3B ), body length (ρ = 0.78, P = 2.5 × 10 −8 ; P PGLS = 1.8 × 10 −9 ; Fig 3C ), or ω (ρ = 0.49, P = 2.9 × 10 −5 ; P PGLS = 1.7 × 10 −7 ; Fig 3D ). Similar results were obtained for the testis and ovary for BUSCO genes ( Fig 2E ) or all protein-coding genes ( S2E Fig ).

Although ATI is ideally assessed by CAGE-seq data [ 58 ], such data are available for only several species in our collection ( S1 Data ). Because ATI can be inferred from RNA-seq data [ 59 , 60 ], we chose to use RNA-seq to compare ATI across species to allow the inclusion of a broader range of species in our analysis. We used SEASTAR, a method known to outperform other tools [ 60 ], to predict ATI sites. To validate the SEASTAR prediction, we collected a set of samples sequenced by both CAGE-seq and RNA-seq [ 61 ] and, respectively, identified ATI sites from CAGE-seq and from RNA-seq using SEASTAR for all protein-coding genes. Analogous to the APA analysis, we quantified the expression level for each ATI site, the total percentage usage of minor ATI sites for each gene, and the average total percentage usage of minor ATI sites per gene in a species (see Materials and methods ), and found them to respectively exhibit a significant, positive correlation between estimates from CAGE-seq and those from RNA-seq (ρ > 0.32, P ≤ 0.02; S2A – S2C Fig ). These findings support the reliability of predicting ATI site usage from RNA-seq data.

Next, we used RNA-seq data to estimate the average total percentage usage of minor APA sites per BUSCO gene in each of 75 species. We started with the brain because the number of species with RNA-seq data from this tissue is the highest. Across species, the above estimate of species transcript diversity reduces with N e (ρ = −0.73, P = 4.3 × 10 −4 ; P PGLS = 1.8 × 10 −2 ; Fig 2B ), but increases with life span (ρ = 0.61, P = 9 × 10 −5 ; P PGLS = 4.2 × 10 −2 ; Fig 2C ), body length (ρ = 0.70, P = 2.4 × 10 −6 ; P PGLS = 4.4 × 10 −7 ; Fig 2D ), and ω (ρ = 0.63, P = 8.6 × 10 −9 ; P PGLS = 3.5 × 10 −3 ; Fig 2E ). Similar results were observed for the ovary and testis ( Fig 2F ). When the above transcript diversity in a species was calculated using all protein-coding genes, the patterns remain qualitatively unchanged ( S1F Fig ). Hence, patterns of interspecific variation in transcript diversity caused by APA support the error hypothesis.

We previously reported a negative correlation between the total percentage usage of minor APA sites of a gene and the gene expression level in each of five mammals studied, supporting the error hypothesis [ 32 ]. To assess whether this pattern extends beyond mammals, we repeated the analysis using APA sites predicted from RNA-seq data and found that the negative correlation persists in 163 of the 166 samples analyzed ( S1E Fig ), suggesting that this pattern is generally true across animals.

To increase the number of species in the APA analysis, we predicted APA sites using RNA-seq data [ 48 , 54 ]. In particular, TAPAS leverages the Pruned Exact Linear Time algorithm, RNA-seq data, and gene structure information to predict APA sites and their abundances [ 55 ] and is known to outperform other tools [ 56 ]. We validated the performance of TAPAS by comparing the APA sites identified by 3′-end-seq with those predicted by TAPAS from RNA-seq using a dataset in which CD4 T cells were sequenced by both 3′-end-seq and RNA-seq [ 57 ]. We quantified APA site usage levels in both data types (see Materials and methods ) and found a significant positive correlation between them (ρ > 0.27, P < 0.01; S1B Fig ). We also computed the total percentage usage of minor APA sites for each gene in both data types and again observed a significant positive correlation between them (ρ > 0.16, P < 0.01; S1C Fig ). To validate the applicability of APA site prediction by TAPAS in multiple species, we acquired RNA-seq datasets from the 10 species with 3′-end-seq data; while 3′-end-seq and RNA-seq were not generated from the same samples, they were from the same tissues in each of these species ( S1 Data ). We found that the average total percentage usage of minor APA sites per gene is significantly correlated between 3′-end-seq and RNA-seq data across species (ρ = 0.79, P = 9.8 × 10⁻³; S1D Fig ). These results confirm the reliability of predicting APA sites from RNA-seq data.

(A) Relationship between N e and transcript diversity caused by APA estimated from 3′-end-seq. (B–E) Relationship between transcript diversity caused by APA estimated from RNA-seq and N e (B) , life span (C) , body length (D) , or ω (E) in the brain. P-values from Spearman’s correlation and PGLS are shown. N represents the number of species included. (F) Correlation between transcript diversity caused by APA and N e (or proxies) in three tissues. BUSCO genes are used in all panels. The data underlying this Figure can be found in https://doi.org/10.5281/zenodo.18514977 .

Accurate detection of APA typically relies on 3′-end RNA sequencing (i.e., 3′-end-seq), which can identify precise APA sites and their relative usages [ 48 ]. Several library preparation methods are available for 3′-end-seq [ 49 ], such as 3′READS [ 50 ], 3P-seq [ 51 ], and PAS-seq [ 52 ]. However, only a limited number of species have 3′-end-seq data [ 53 ]. We collected 3′-end-seq datasets from 10 species with N e ( S1 Data ). For a given gene, let us refer to the most frequently used APA site as its major APA site, which is likely to be functionally the best APA site, and all other APA sites as its minor APA sites. We measured the transcript diversity of the gene caused by APA by the total percentage usage of its minor APA sites, which equals the total 3′-end reads of its minor APA sites divided by the total APA reads of the gene. We then averaged the transcript diversity across all genes considered in a species to represent the overall transcript diversity due to APA for the species. We found a negative correlation between transcript diversity and N e across species, regardless of whether we considered all protein-coding genes ( S1A Fig ) or only BUSCO genes ( Fig 2A ). However, the negative correlations (or those from PGLS) were not significant, potentially due to the limited number of species in the analyses.

Of the 100 species, 26 have published N e ( Fig 1C and S2 Data ). Given that most species lack published N e , we resorted to three other parameters as N e proxies: the body length, life span, and nonsynonymous to synonymous substitution rate ratio (ω). Because larger animals and longer-lived animals tend to have smaller N e , body length and life span have been used as N e proxies [ 36 , 44 , 45 ]. Under the nearly neutral theory and the neutral assumption of synonymous mutations, ω is expected to decline with N e [ 41 ] so can also be a proxy for N e . However, when synonymous mutations are frequently non-neutral as has been documented in some species [ 46 ], the validity of the above expectation is uncertain; we therefore examined it empirically (see below). We collected body length and life span data from previous studies [ 36 , 47 ] and estimated ω for each species (see Materials and methods and S3 Data ). We correlated the three proxies with N e and confirmed their negative correlations ( Fig 1D – 1G ), as reported previously [ 36 ]. Therefore, the error hypothesis predicts that transcript diversity should decline with N e or increase with the three N e proxies considered here. Note that throughout this study, we used two types of correlation analysis. The first is the simple rank correlation, whereas the second is Phylogenetic Generalized Least Squares (PGLS) regression, which controls for the phylogenetic relationships of the species considered in the data (see Materials and methods ).

Because transcript diversity estimated from different tissues may not be comparable across species, we focused on three tissues (brain, ovary, and testis) best represented in the transcriptomic data, respectively covering 67, 51, and 48 species ( Fig 1B and S1 Data ). In interspecific comparisons, we initially analyzed all protein-coding genes in each species ( Fig 1 and S1 Data ). However, gene set variations among species could introduce a confounder in our comparison. We therefore focused on Benchmarking Universal Single-Copy Orthologs (BUSCO) genes in an additional comparison ( Fig 1B and S1 Data ), as was done in previous studies [ 36 , 43 ].

(A) Phylogenetic tree of the 100 metazoan species considered. (B) Available genome and transcriptome datasets of each species concerned, including CAGE-seq, 3′-end-seq, RNA-seq, coding genes, and BUSCO genes. (C) N e and proxies. Spearman’s correlation and Phylogenetic Generalized Least Squares (PGLS) regression between N e and life span (D) , body length (E) , and the nonsynonymous to synonymous substitution rate ratio ω computed using all BUSCO genes (F) and all one-to-one orthologous genes (G) across species. Each dot represents a species, colored according to its clade in (A) . The data underlying this Figure can be found in https://doi.org/10.5281/zenodo.18514977 .

Based on two recent studies [ 36 , 43 ], we assembled a list of 100 diverse metazoan species with available genome sequences ( Fig 1A ). We collected publicly available genome, genomic annotation, and transcriptome data from these 100 species ( Fig 1B and S1 Data ). The transcriptomic data are the basis of transcript diversity estimation and they comprise Cap Analysis of Gene Expression Sequencing (CAGE-seq) data from seven species for direct measurement of ATI, 3′-end-seq data from 10 species for direct measurement of APA, and RNA-seq data from 75 species for detection of AS and prediction of ATI and APA ( Fig 1B and S1 Data ).

Discussion

Transcript diversity primarily arises from ATI, APA, and AS. Although past studies have provided substantial genomic evidence for the error hypothesis of ATI [33], APA [32,65], and AS [34], these studies focused on a small number of species. As a result of this limitation and an increasing number of reports of cases of functional ATI, APA, or AS, the general biological significance of transcript diversity remains controversial. In the present study, we expanded the analysis to 75 species and showed that the previous finding from a small number of species generally hold across animals. More importantly, we found that the transcript diversity of a species declines with the species’ N e or its proxies, as predicted by the error hypothesis.

A central theoretical underpinning of the non-adaptive paradigm relevant to our results is the drift-barrier model, which was first proposed by Lynch [66] to explain the mutation rate variation across species. This model predicts that the efficacy of selection is limited due to genetic drift and mutation bias, such that phenotypic traits of species with smaller N e , where drift is more potent, are less optimized than those of species with larger N e [42]. For transcript diversity, the drift-barrier model predicts that errors in transcriptional and post-transcriptional processing (e.g., incorrect AS, imprecise polyadenylation, or aberrant transcription initiation) that generate non-functional or weakly deleterious transcript isoforms will persist at higher frequencies in species with smaller N e . Our observation that transcript diversity declines with N e (and its proxies) confirms this prediction and hence supports the drift-barrier model.

Comparative studies across species often encounter confounding factors that could bias the outcome. For instance, estimating transcript diversity in this study relies on transcript annotations, which vary in completeness across species, with model or well-studied organisms typically having more comprehensive annotations, potentially introducing an interspecific bias. To minimize this bias, we performed de novo transcript assembly for each species using RNA-seq data, ensuring uniform annotation processes and quality (see Materials and methods). Although our approach may reduce annotation quality for model organisms, it mitigates the potential interspecific bias. Indeed, the observed patterns are unaltered by including (Figs 2E, 3D, and 4D) or excluding (S4 Fig) model species. Similarly, genome size and complexity can influence genome annotations and transcript diversity assessments. To address this issue, we employed BUSCO genes. These genes are single-copy highly conserved orthologs that are unaffected by genome size or complexity across species, ensuring comparability among taxa. Indeed, when using BUSCO genes to compute transcript diversity, we found that our conclusions hold regardless of whether genome size and complexity are controlled or not (S5 Fig). Thus, our results are robust to the above potential confounding factors in multispecies comparisons.

As mentioned, Benitiere and colleagues (2024) also reported a negative correlation between N e and transcript diversity caused by AS across 53 species [36]. Nevertheless, our methodology differs from Benitiere and colleagues’s in several aspects. First, we calculated transcript diversity at the gene level, whereas Benitiere and colleagues measured it at the intron level. Second, Benitiere and colleagues combined multiple RNA-seq datasets to detect splicing events, which improved splicing event detection but introduced heterogeneity among tissues. By contrast, we compared the same tissue across species, which made the interspecific comparison fairer. Third, Benitiere and colleagues analyzed 53 species, most being insects, while our analysis encompassed 75 animals with a broader phylogenetic sampling. Fourth, in addition to the three N e proxies used by Benitiere and colleagues in the correlation analysis, our study also used N e from 26 species. Notwithstanding these methodological differences, the findings of the two studies are consistent.

In the debate about the biological significance of AS, several authors noted a significant positive correlation between the amount of AS of a species and its organismal complexity measured by the number of cell types [67,68]. Chen and colleagues [68] reported that this correlation remains even after the control for the species’ N e , suggesting that AS is adaptive and is at least partially responsible for organismal complexity. Their analysis was recently criticized by Benitiere and colleagues [36] for using nucleotide diversity at synonymous sites (π S ) as a proxy for N e . Synonymous mutations are often non-neutral [46], but even when they are neutral, π S is determined by both N e and the mutation rate per site per generation, the latter of which varies across species [42]. Hence, π S is not an appropriate proxy for N e . While Benitiere and colleagues suggested that ω would be a more appropriate proxy for N e , they did not perform the actual partial correlation analysis. We therefore investigated the relationships among transcript diversity caused by AS, number of cell types, and ω across 12 species for which all three estimates are available in our brain tissue dataset (S4 Data). Consistent with the finding of Chen and colleagues [68], we observed a positive correlation between the transcript diversity caused by AS and the number of cell types across species even after the control for ω, but this partial correlation did not reach statistical significance in the PGLS analysis (S6 Fig). That is, after the control for the phylogenetic relationships in the data and ω (as a proxy for N e ), there is no significant partial correlation between organismal complexity measured by the number of cell types and transcript diversity caused by AS. We note that, even if the above partial correlation is significant, it does not mean that AS underlies organismal complexity. This is because, to demonstrate that AS contributes to organismal complexity, one needs to show that AS varies among cell types and plays a role in the functional diversity among cell types, which will be an interesting direction to pursue in the future when AS can be reliably assessed from single-cell RNA-seq data.

It is important to note that beyond APA, ATI, and AS, there are other variations in gene expression that result in gene product diversity, including post-transcriptional modifications (e.g., RNA editing [69] and m5C modifications [37]), translation variations (e.g., alternative translation initiation [70], mistranslation [71], and stop-codon read-through [72]), and post-translational modifications (e.g., phosphorylation). Notably, a recent study of the mis-transcription rate reveals a narrow range of variation across the tree of life [73]. Hence, it would be highly valuable to conduct cross-species comparisons as performed here for additional types of transcript diversity when appropriate data become available from a sufficient number of species.

[END]
---
[1] Url: https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3003671

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/