Evolution of Exon-Intron Structure and Alternative Splicing

Despite significant advances in high-throughput DNA sequencing, many important species remain understudied at the genome level. In this study we addressed a question of what can be predicted about the genome-wide characteristics of less studied species, based on the genomic data from completely sequenced species. Using NCBI databases we performed a comparative genome-wide analysis of such characteristics as alternative splicing, number of genes, gene products and exons in 36 completely sequenced model species. We created statistical regression models to fit these data and applied them to loblolly pine (Pinus taeda L.), an example of an important species whose genome has not been completely sequenced yet. Using these models, the genome-wide characteristics, such as total number of genes and exons, can be roughly predicted based on parameters estimated from available limited genomic data, e.g. exon length and exon/gene ratio.


Introduction
Recent advances in high-throughput DNA sequencing led to significant progress in complete genome sequencing and opened unprecedented opportunities for comparative genome studies [1][2][3][4]. The complete genome sequences are publicly available from constantly growing databases, such as the National Center for Biotechnology Information (NCBI) GenBank, and can be readily analyzed and compared for a number of evolutionarily distant species. The early comparisons revealed that the number of genes and metabolomic complexity progressively increase as species become more evolutionarily advanced [5][6][7][8], but their anatomical, morphological, physiological and behavioral complexity does not linearly correlate with the total number of genes discovered. For instance, whereas the number of protein coding genes in the human genome is only 14% greater than in the roundworm Caenorhabditis elegans, the evolutionary differences between these two species are immense. This suggests that regulatory and posttranscriptional processes might play an increasingly more important role throughout evolution. There are numerous mechanisms, processes and structures that affect gene regulation, such as methylation, chromatin structure, regulatory elements, transcription factors, polyadenylation, post-translational modifications and compartmentalization of proteins, and others (for review see [9]). RNA editing and alternative splicing (AS) are two distinct posttranscriptional processes that can, however, increase proteomic complexity and number of various proteins without increasing the number of genes. While RNA editing, particularly common in organelles (for review see [10,11]), can lead to nucleotide insertion or deletion, in higher Eukaryotes base modifications are prevailing [11][12][13]. In the process of AS additional variants are created among the mature mRNA transcripts through modification and rearrangement of exons (e.g. [14,15]). Both mechanisms promote adaptive and evolutionary potential of species without increasing the number of genes and maintenance cost that could be associated with it. For instance, due to AS the total hypothetical number of various proteins encoded by the DSCAM gene can reach 38,016 in Drosophila melanogaster [16]. Therefore, one may expect that more evolutionarily advanced organisms have more elaborated and complex AS. We addressed this hypothesis in more detail in our study. Our objectives were to examine exon-intron structure in genomes of completely sequenced and fully annotated species, to infer AS data and to use this information for defining relationships between genes and proteomic complexity. We expect that these relationships can be used to predict the anticipated exon-intron structure and proteomic complexity in non-model species with large genomes, such as pines, that may remain unsequenced for a while. We applied our findings to loblolly pine (Pinus taeda L.), one of the most-studied coniferous species, which has a very large genome of 24.56 pg (,24 Gb) [17]; complete genome sequencing for loblolly pine is underway [18], but is still problematic and unavailable. The obtained knowledge is also essential for understanding the genetic control of the metabolomic complexity and functionality in the studied species and the evolutionary significance of AS in general.

Analysis of complete genomes
First, to explore general trends in the completely sequenced genomes, we analyzed basic statistics (Table 1 and Table S1). Previously Lynch and Conery [6] showed strong positive correlation between genome size and gene number in multiple species. Here, we observed an increase in the number of genes (Fig. 1A), gene products and total number of exons as species advance evolutionarily. The exon/gene ratio also increases (Fig. 1B), but the mean and median exon length becomes shorter (Fig. 1C), whereas CDS length remains relatively constant. Because the changes in parameter values demonstrated clear trends that followed evolutionary advancement, we proceeded with an in-depth correlation analysis. The results of regression analysis and estimates of the parameters are summarized in Table 2 and Table S2.

Average exon length as a predictor
Generally, the length of exon can be approximated from a limited sample of genes. Hence, we explored its potential for predicting other genomic parameters in incompletely sequenced species. A very strong negative correlation was observed between mean exon length and total number of exons (r 2 = 0.937, r 2 adj = 0.936; Fig. 2A). The negative correlation was weaker but statistically significant between mean exon length and either number of protein coding genes (r 2 = 0.712, r 2 adj = 0.703; Fig. 2B) or the total number of genes (r 2 = 0.706, r 2 adj = 0.697, Fig. S1). Mean exon length also strongly negatively correlated with exon/ gene ratio (r 2 = 0.957, r 2 adj = 0.956; Fig. 2C). In all these cases the correlations were not linear. No statistically significant correlation was observed between mean exon length and mean CDS length (P = 0.132; Fig. S2). These results proved that exon length is, indeed, a robust predictor.
Similar correlations were observed for median exon length, and the r 2 values were close to those obtained for mean exon length (see Table 2 for details).

Number of genes and exons in genomes
Since both number of genes and number of exons increase with the species complexity, we measured how robust the correlation between these parameters was. Strong and nonlinear but positive correlations were found between the number of genes and number of exons (r 2 = 0.897, r 2 adj = 0.894 for protein coding genes, Fig. 2D; r 2 = 0.891, r 2 adj = 0.888 for the total number of genes, Fig. S3). Also very strong positive correlation was observed between the total number of genes and number of protein coding genes (r 2 = 0.996, r 2 adj = 0.996; Fig. 2E). These tight relationships show clear directions in the genome evolution, where not only amount of genetic information is increasing but also is accompanied by fragmentation that facilitates its dynamic use. The linearity and extremely high r 2 value between all genes and protein coding genes also shows that the non-protein coding gene fraction changes proportionally.

Exon/gene ratio as a predictor
Similarly as the exon length, exon/gene ratio has a high predictive power for complex genomic parameters. Strong positive linear correlation was observed between exon/gene ratio and the total number of exons in the genome (r 2 = 0.864, r 2 adj = 0.860; Fig. 2F). The relationship between exon/gene ratio and the number of protein coding genes was also strongly positive but rather nonlinear (r 2 = 0.648, r 2 adj = 0.638; Fig. S4).

Alternative splicing
The traditional methods of using EST data to study AS are sensitive to the EST coverage [19]. To diminish biases related to this problem we used the NCBI GenBank annotations of the wellannotated genomic data for 12 model species in our study. Despite numerous AS studies (see Discussion) we are unaware of any that would take an advantage of thorough annotations. The five most annotated species (A. thaliana, C. elegans, D. melanogaster, M. musculus and H. sapiens) were analyzed in more detail (Table 3). Among the five most common AS types, alternative 39 splice sites type (A3) was the most frequent in A. thaliana and C. elegans, but exon skipping (ES) was the most frequent in D. melanogaster, M. musculus and H. sapiens.
The frequency of the ES type increases following organism complexity and reaches more than 51% of all AS types in human (Fig. 3). It is accompanied by a decrease in the intron retention (IR), A3 and alternative 59 splice sites (A5) types. In the plant model, A. thaliana, the most common type is A3 (44.3%), followed by IR (19.7%).

Alternative splicing ratio
The AS ratio, that is the ratio of the alternatively spliced genes and the total number of protein coding genes was highest (0.186) in Pan troglodytes among all analyzed species. When only the mostannotated species were considered, the highest AS ratio was observed in human (0.080), followed by the one in D. melanogaster (0.073) and A. thaliana (0.066; Table 3). In general, the ratio increases with evolutionary progress (Fig. 1D). Although these estimates seem low as compared to other AS studies (see Discussion for details), we may expect a rapid growth in the data abundance in the near future due to the fast development of novel genomic techniques (e.g. [20]), which will help fully understand the effects of this mechanism.
AS negatively correlated with exon length and occurred more frequently in organisms with shorter exons ( Among the five most-annotated species, both Shannon's diversity index and equitability were highest in A. thaliana (H = 7.38, E = 0.985; Table 3) showing high richness and evenness of distribution. Second high value (H = 7.26) was observed in human, but the evenness was lower (E = 0.966).

Predictions for other species with large genomes such as Pinus taeda
Having defined regression equations between the key genomic parameters, we applied these models to P. taeda, an important but incompletely sequenced species. Based on the 99 complete CDSs available in GenBank, the mean and median transcript lengths were practically the same (1278 bp; Table 4). However, the mean and median exon lengths were very different -334.8 and 198 bp, respectively. These estimates are very preliminary and based only on 21 exons. An additional 43 complete exons were identified in partial CDSs. Their length was shorter -166.2 bp on average, but these estimates could be biased toward shorter exons due to the PCR-biased amplicon resequencing. It is worth noting that any Table 2. Predicted values for exon-intron gene structure and alternative splicing (AS) parameters for an organism with mean and median exon lengths of 334. 8 and 198.0 bp, respectively, such as observed in Pinus taeda, based on results of regression analysis.

Response (y)
Factor (x) R 2 R 2 adj P-value at 95% CI Figure  species, completely or incompletely sequenced, with a known estimate of the exon length could be investigated here in place of P. taeda. Moreover, the accuracy of the exon length estimation for any species will increase if more data are available and if the ascertainment bias (e.g. resulting from overrepresentation of certain gene families in the sample, or from inclusion of only partially sequenced genes) is kept at the minimum. Based on the regression models created for the 36 complete genomes, we computed estimates for a hypothetical species with an average exon length such as the one observed in P. taeda ( The estimates differed slightly when median exon length was used (see Table 2 for details).

Discussion
Despite significant progress in sequencing technologies, complete genomic data are still limited for eukaryotic organisms; and, more importantly, only a few extensively studied model species have been well annotated and featured. The most abundant data have been collected for microbial, fungal and some animal genomes, while vascular plants have been understudied. Only a few species have been completely sequenced and annotated in this underrepresented group so far, such as Arabidopsis thaliana, Oryza sativa, Vitis vinifera, Physcomitrella patens, and, recently, Populus trichocarpa, which remains relatively poorly annotated and featured. Therefore, due to insufficient experimental data, it is very likely that not all gene transcripts and AS products have been recorded in GenBank even for the best-studied species; this can cause underestimation of AS ratios in our study. As more experimental data are collected, the situation gradually improves with every new genome build that updates the number of genes, exons, and their locations on the chromosome. AS has different types, occurs at different developmental stages and tissues, and can be affected by environmental factors [21,22]. AS is still insufficiently studied, and, therefore, not all AS events and types are well documented in the databases. Moreover, precise inference is more difficult due to   incomplete annotation, along with different stringency criteria, customary thresholds to classify true and erroneous AS events, and various gene models used in different species. To avoid these complications, we limited our analyses only to extensively studied completely sequenced genomes. However, we hope that the results obtained can be used for predictions in insufficiently studied and incompletely sequenced organisms, such as pine in our study, in which investigation of exon-intron architecture is impeded by the availability of the data.

Alternative splicing ratio
In general, both AS and number of genes are higher in more evolved organisms (Figs. 1A, D). However, surprisingly the AS rate was not always as high in more evolutionarily advanced species as expected and did not correlate linearly with their evolutionary progress. For instance, the ratio observed in D. melanogaster (0.073) was higher than in A. thaliana (0.066), but very close to the one in human (0.080). This could suggest that a relatively small number of genes in D. melanogaster compared to human (14,807 vs. 25,074, respectively) is compensated by a higher AS rate that increases proteomic and metabolomic complexity. We cannot completely exclude that the discrepancies are explained by insufficient AS data, but we can conclude in general that both the number of genes and the AS rate increase in more evolutionarily advanced species. Certainly, more experimental data are needed to increase the precision of estimates and predictions of AS ratios. For instance, a number of studies demonstrated AS in rice [23][24][25], but it is not documented in rice genomic data from the NCBI GenBank database, a likely indication that annotation of the rice genome is still in progress.
A substantially higher AS rate for shorter exons found in the study (Fig. S5) is consistent with previous studies that suggested that both tandem exon duplication [26] and insertion of noncoding intron sequences [27] could promote AS. Both within-gene duplication and AS would have less drastic effect on functionality of final proteins when they both deal with mutually exclusive exons (ME) that have shorter lengths. In addition, converting a part of an intron into an exon via AS has the risk of including a stop codon. This risk is higher when alternative exon sequences are longer.
Although not as drastic, higher exon/gene ratio is also associated with a higher AS rate (Fig. S6). This can be observed from the above described correlation of exon length and AS because more exons per gene mean both shorter exons (assuming a constraint on the final gene product length) and more options for AS. Exon/gene ratio, similarly to AS, increases in evolutionarily advanced species (Figs. 1B, D). More advanced species also show higher numbers of CDSs. AS increases as number of CDSs increases (Fig. S7), playing an important role in creating higher proteomic complexity.
Previous studies reported the ES type of AS as the most frequent in mammals [28,29]. Our results are consistent with these findings. ES accounted for 51.0% of AS in human and 41.3% in mouse. Moreover, frequency of ES increases with complexity. Kim et al. [30] used a modified approach that required the final number of ESTs in the compared organisms' genes to be the same to mitigate the bias in data availability for the studied species. They also found a high frequency of ES in mammals (,40%) and low IR (,10%).
It is the opposite in plants, where IR type was the most frequent (over 50%) in both rice and thale cress [24,31]. The A3 was the second most frequent type, while ME was the least frequent type. Nagasaki et al. [28] reported that IR accounted for over 42% of AS events in thale cress and 55% in rice. Although our analysis also found a very low level of ME in thale cress (0.2%), the most abundant type was A3 (44.3%) followed by IR (19.7%). Assuming that none of the classes is underrepresented in the dataset we used, this could indicate that IR tends to be overestimated in the EST/ cDNA based studies, possibly due to the highest incidence of nonsense-mediated mRNA decay (NMD)-targeted products in this class. Wang and Brendel [31] estimated that ,43% of AS events in Arabidopsis are potential NMD candidates, with IR showing the highest incidence of 40-48%. Conversely, Kim et al. [30] found the rate of IR in Arabidopsis (,30%) less than A3 (,40%). In their study, ES accounted for approximately 5% of all types. McGuire et al. [25] found that in A. thaliana IR accounted for 38.7% of splice variants, only slightly more than A3 (36.8%) and ES (7.7%) events. These results show great sensitivity to the methods and assumptions used. McGuire et al. [25] discussed how including unspliced alignments may affect the outcomes.

Alternative evolutionary scenarios for plants
This study demonstrates that plants and animals may have used different mechanisms and strategies for developing proteomic and metabolomic complexity. The AS rate is low in plants compared to animals, whereas the number of genes is high (Figs. 1A, D). This could indicate that animals have evolved a more efficient system of managing the genomic information that allows them to increase proteomic complexity with the same or smaller number of genes. Flowering plants could have relied primarily on duplications (from exon shifting to entire chromosome or genome duplications that are common in flowering plants) [32], duplication modifications and divergence. In contrast, large genomes in genus Pinus (class Coniferopsida) might be a result of retrotransposon expansion rather than polyploidy [33]. Cui et al. [32] found no evidence of recent genome duplications in P. taeda nor P. pinaster, and Grotkopp et al. [17] estimated that the genome sizes varied from 22.10 pg to 36.89 pg in pines, with the putative common ancestor's genome of 32.09 pg. Nevertheless, sporadic polyploidy has been observed in gymnosperms (for review see [34]).
In order to check how well our regression models fit the real data, we compared the values observed in the three plants in our dataset with those predicted by the models. The total numbers of exons predicted based on the observed average exon length in Oryza sativa and Arabidopsis thaliana were underestimated in both cases (predicted 92,581 and 102,910 vs. observed 128,267 and 138,876, respectively), but the observed values were only slightly greater than the 95% confidence interval upper limit at the individual level (126,092 and 136,543, respectively; Fig. 2A). The observed number of exons in the primitive alga Ostreococcus 'lucimarinus' (9,767) was close to the predicted number (10,021) and fell within 95% CI.
Similarly, the observed numbers of protein coding genes (Fig. 2B) in both higher plants (26,977 in thale cress and 26,777 in rice) were only slightly higher than the upper 95% CI limit at individual level (26,282 and 25,445, respectively), and predicted values were underestimated (17, 694 and 16,888, respectively). The total observed number of genes (Fig. S1) for the two species, again, fell only slightly above the 95% individual level CI (28,245.  27,591 for thale cress and 29,102.26,715 for rice; predicted  values: 18,480 and 17,637, respectively). In both cases the numbers observed in the alga fell within the 95% CI.
The discrepancy between Viridiplantae and other kingdoms can also be observed in the relationship between exon count and number of genes (Fig. 2D, Fig. S3), as well as exon/gene ratio and number of protein coding genes (Fig. S4). In our models the values for the two higher plants fall outside the upper 95% CI at the individual level in all these cases. Interestingly, comparison of the two evolutionarily youngest groups in kingdoms Metazoa and Viridiplantae reveals much longer exons in plants (236.8 bp in A. thaliana and 250.2 bp in O. sativa) than in mammals (ranging from 169.2 bp in Canis lupus familiaris to 179.4 bp in M. musculus), demonstrating that the processes that reduce exon length have been slower in plants. Body plan complexity is much greater in mammals than in angiosperms, and shorter mammalian exons coupled with lower number of genes could indicate greater pressure towards efficient use of gene space. The highest values of Shannon's index and equitability observed in A. thaliana (Table 3) indicated more even AS distribution than in four animal species, despite their higher evolutionary position. Perhaps AS is not the main mechanism in achieving the observed complexity level in plants. If AS is correlated with exon length, then longer exons in plants can imply that the pressure for greater AS rates is not as strong as in the case of higher animals. Other mechanisms, such as more frequent duplications and elevated retrotransposon activity in plants could be responsible for the high number of genes, and also greater exon lengths, through intron loss. Indeed, studies on animal species have shown a negative correlation between gene family size and AS frequency [35][36][37]. This could explain not only lower rates of AS observed in plants but also the different patterns of AS forms, potentially increasing chances of NMD, a phenomenon not very well studied in plants [23]. Conversely, building upon the theory proposed by Lynch and Conery [6], Babenko et al. [38] suggested that intron gain/loss is not a commonly ongoing process, but rather may be triggered by certain dramatic evolutionary events that lead to long-term bottlenecks. Therefore the observed differences in exon lengths could be merely due to chance of the ancestors being affected by drastic events in the past. These conclusions seem to be supported by Sammeth et al. [29], showing rather abrupt differences between invertebrates and vertebrates.
Current genomic data are insufficient to build separate robust regression models for plants. The conclusions about the total number of exons, the total number of genes and the total number of protein coding genes in P. taeda may therefore be biased; and the true values may be close to the upper 95% CI limit on the individual level, higher than the ones predicted by the proposed models (see below).

Short exons promote genomic complexity
The strong relationship between exon length and total number of exons ( Fig. 2A) as well as exon length and exon/gene ratio (Fig. 2C) suggest that shorter exons increase potential for AS. Indeed, a much higher ratio of AS was observed in organisms with shorter exons (Fig. S5). However, more evolutionarily advanced organisms not only have shorter exons, but also more genes (Fig. 2B, Figs. 1A, C, and Fig. S1). The presence of shorter exons increases the potential for exon shuffling along with exon duplications; and, as a complement of AS, both increase proteomic and metabolomic complexity. It is likely that both evolved simultaneously and synergistically to amplify their effects on increasing physiological, behavioral and morphological complexity of the organisms through positive feedback loop-like mechanisms.
No statistically significant correlation was found between exon length and CDS length (Fig. S2). Since 3-dimensional protein structure and binding sites determine protein functionality, the length of the coding sequence seems to be of primary importance, and therefore the variation in the transcript length may be constrained. This could suggest that in the process of evolution, partitioning of the ancestral coding sequences has been occurring rather than extension through e.g. hypothetical stacking of coding blocks together. Such a process could have stimulated splicing out duplicated exons, eventually leading to alternatively spliced forms.
At the genome level, most of the species with less than 10,000 genes had a very small number of exons ( Fig. 2D and Fig. S3). Consequently, the number of exons per gene was low in these species (Fig. 2F, Figs. 1A, B and Fig. S4). These observations show a general trend of genomic complexity increasing in evolutionarily advanced species.
The shortest exons identified in some of the analyzed species (including three plants) were only 1 bp long. We did not find any peer-reviewed publications experimentally confirming this observation. In previous studies Long et al. [39] identified single base pair exons, and Deutsch and Long [40] identified exons as short as 1 amino acid in a number of species including A. thaliana and human, although the exact length in bp is not clear. An experimental approach is necessary to find support for these structures and to verify that it is not an artifact resulting from exon/intron model assumptions. For instance, Kondrashov and Koonin [26] used 9 bp as threshold.

Implications for Pinus taeda
Due to the small sample size, the P. taeda exon length estimate may be significantly biased. An alternative would be to include in the study completely sequenced exons from only partially sequenced genes. However, in this scenario shorter exons would be overrepresented due to the PCR amplicon length bias (typically a few hundred bp), making the mean and median underestimated. The observed exon/gene ratio was 4.000 based on 5 CDSs and 20 exons. This estimate is very close to the predicted exon/gene ratio of 4.245, based on the regression model, when the average exon length was the predictor or 3.658 when the median exon length was used ( Table 2). The predicted values based on the average exon length for the other two higher plants analyzed were also close to the observed values (5.970 vs. 5.148 in thale cress and 5.654 vs. 4.790 in rice). The observed values for all three species fell inside the 95% CI at the individual level.
The number of protein coding genes and total number of genes expected in an organism with the average exon length of 334.8 bp (such as in P. taeda) is 13,871 and 13,288, respectively. These values seem to be underestimated as far as P. taeda is concerned, especially when compared with the other analyzed plants. Moreover, there are currently about nineteen thousand unique sequences in the NCBI UniGene database for P. taeda. The model severely underestimates the number of genes in the other two vascular plants described above as well. The number of protein coding genes in A. thaliana is underestimated by about 34.4% and O. sativa by about 36.9%. If P. taeda followed this bias, and the expected number of protein coding genes was also underestimated by approximately 35%; that would mean about 7,155 underestimated genes, which would raise the predicted number of protein coding genes to about 20,443 in this species, making this number more realistic.
Regression models that are based exclusively on the three examined plants and that follow the same logic as in the case of the 36 studied genomes also demonstrated that higher numbers are expected for loblolly pine ( Table 2). The total number of genes expected would be 20,923 (upper limit is 45,975 at 95% confidence level; model significant at 93.7% confidence level) and the number of protein coding genes 19,785 (upper limit is 32,508 at 95% confidence level). These numbers seem to be more realistic when compared to the observed values in other higher plants, especially considering broad confidence intervals. Although the above findings for loblolly pine are restricted by the limited availability of the data, this study is the first attempt to address the questions of exon-intron structure and genomic complexity in this species, and will likely stimulate further studies.

Conclusions
This study confirmed the general trend of increasing number of genes, gene products, and exons in the genome, along with higher exon/gene ratio and AS ratio as species become more evolutionarily advanced. We demonstrated that parameters easily computable from small data samples (e.g. exon length or exon/gene ratio) are relatively good predictors of characteristics that are difficult to assess, such as total number of genes, gene products and exons. We also showed that taxonomic kingdoms may require different model calibration as their strategies to increase complexity throughout evolution have been different. As more genomic data become available and more species representing various taxonomic groups are annotated, these models can be tuned or applied to specific monophyletic groups, which will improve precision of the predictions.

Selection of completely sequenced species for analysis
We selected the 36 most-annotated and featured species (Table 1) from eukaryotic genomic assemblies available in the NCBI GenBank [41]. For 10 species in our set (Arabidopsis thaliana, Caenorhabditis elegans, Drosophila melanogaster, Encephalitozoon cuniculi, Eremothecium gossypii, Homo sapiens, Mus musculus, Oryza sativa, Saccharomyces cerevisiae and Schizosaccharomyces pombe), the annotation involved at least partial curation of the records. AS in rice was not documented in the present NCBI GenBank genome annotation, although it has been reported previously [23,24]. Therefore, this species was excluded from the AS analysis.

Source of data
All genomic data were downloaded from the FTP directory of the NCBI GenBank (ftp://ftp.ncbi.nih.gov/genomes/MapView/). Sequences for P. taeda were downloaded from the Nucleotide database from the NCBI GenBank.

Genomic data analysis
Genomic data were analyzed using Perl scripts specifically written for this study (available at http://treenome.tamu.edu/). The downloaded files were screened, and chromosome ID, position and orientation of the exons on the chromosome, feature ID, AS type, transcript accession number, and group label were traced, partitioned and analyzed. Pseudogenes, mitochondrial, plastid and insufficiently annotated genes were excluded from further analysis. Total numbers of genes, protein coding genes and their coding sequences (CDSs) were calculated for each species. The number of exons and their boundaries were determined based on the coding structure of each protein coding gene ID recorded in their corresponding CDS section. For each gene supported by more than one CDS, the alternative coding sequences were compared with each other. Cases when corresponding exons had different boundaries or no matching counterpart were qualified as AS variants. Average and median lengths were calculated for both exons and CDSs. The exon estimates were computed based on all unique exons found in the genome. All CDSs, including alternatively spliced forms, were considered for estimation of the average and median CDS lengths. The exon/gene ratio was defined as the average number of exons per protein coding gene. AS ratio was defined as the ratio between the number of alternatively spliced and all protein coding genes. AS variants were categorized using the binary approach described by Nagasaki et al. [42]. Shannon's index H and equitability E were calculated for genes with AS to reflect richness and distribution evenness of AS forms [43].

Parameter estimation for Pinus taeda
In total, 99 complete CDS sequences representing protein coding genes in P. taeda were downloaded from NCBI GenBank, and their CDS structure and length were analyzed. The data were prescreened by a Perl script and rearranged manually. The majority of these sequences represented mRNA/cDNA. Only five CDSs represented genomic sequences and could provide complete information about exon-intron structure. Average and median CDS and exon lengths were calculated based on this information. Using the mean exon length and regression models developed based on the genomic data for other species, we computed the expected total number of exons, total number of genes, exon/gene ratio, and total number of protein coding genes for P. taeda. Software package JMP version 5 was used for the statistical analysis.