Improved structural annotation of the P. patens genome

The cosmoss.org Physcomitrella patens V1.6 genome annotation reported here is the result of iterative rounds of evidence mapping, repeat masking, gene structure prediction, filtering and model selection and harbors annotation of protein-coding genes, transposable and repetitive elements and, for the first time, definition of non-protein-coding loci. The release comprises 32,275 protein-coding genes, 432 tRNA loci, 798 rDNA regions, 229 miRNA precursors (108 families) [51], 213 snRNA genes, and 6 SRP (signal recognition particle) loci. Considering the number of miRNA families, P. patens with 108 families has an intermediate position between the green alga C. reinhardtii (47 families) and the flowering plant A. thaliana (187 families) [51]. Consistent with previous findings [3, 10] about half of the genome consists of full length LTR retrotransposons and related fragments including chromodomain-containing gypsy LTR retrotransposons (Tcn1) shared by fungi and non-flowering plants [56].

V1.6 protein-coding gene predictions are based on multiple sources of evidence. Of prime importance are ESTs (Additional file 1: Table A1) from 19 different experimental conditions, tissue types and developmental stages providing a reliable basis for gene structure prediction. The combined transcript evidence was used to train species-specific prediction models using SpliceMachine [57] and EuGène [58], followed by the generation of weighted consensus gene structures using EVidenceModeler [59] and PASA [60]. The prediction procedure was repeated iteratively. Each round involved several filtering steps for separation of non-protein coding, repeat and transposable element-associated genes. The interim versions V1.3-1.5 were not published but are the basis for V1.6 (see Additional file 2: Table A2). The whole protein-coding gene prediction and annotation process is summarized in Additional file 3: Figure A1. The annotation release V1.6 comprises 32,275 loci coding for 38,357 protein-coding transcripts (see annotation releases overview Table 1). 26,722 (~83%) loci are supported by transcript evidence (i.e. EST or full-length cDNA). The average V1.6 gene has a mean length of 2369 bp and a transcript length of 1389 bp. The number of only 1582 unchanged gene models from V1.1 to V1.6 is an excellent indicator of the extent of changes and improvements that led to the current release. While the changes in release V1.2 (Table 1; [3]) were restricted to the removal/filtering of non-protein-coding genes (7972 filtered models), the complete annotation process leading to V1.6 (Additional file 3: Figure A1) resulted in 22,307 (~80% of V1.2) updated models, 8387 (25% of V1.6) new loci and 4077 (15% of V1.2) models which have been removed due to non-protein coding/transposable element origin. Of the new loci, 1535 transcripts (1338 genes) are part of a gene family with at least one additional plant species, and 2196 transcripts (1456 genes) encode at least one Pfam domain. Also, published genes and gene models released by the scientific community were mapped, manually curated and integrated into V1.6.

Inspection of UTR annotations indicates that gene structure completeness is much improved in V1.6. The number of protein-coding loci with both 5’ and 3’ UTRs were increased from 4515 to 15,757 transcripts. The median transcript length increased from 987 (V1.2) to 1248 bases (V1.6) which can be explained by a higher percentage of annotated UTR regions in this release. V1.6 gene models only contain UTR annotations if they are fully supported by transcript evidence.

To further assess model quality and completeness we compared predicted protein sequences from V1.2 and V1.6 based on the coverage of their respective closest homolog in A. thaliana (Figure 1). The direct comparison demonstrates V1.6 models generally as more complete than those in V1.2: 34.5% of the V1.6 proteins with a corresponding V1.2 model cover a higher proportion of the length of their best A. thaliana BLASTP hit.

Annotation and characterization of alternative splicing

We used the PASA pipeline [60] to study the extent and characteristics of alternative splicing (AS) in the moss. We used GenomeThreader [61] transcript alignments in exchange of PASA’s standard GMAP alignments [62] as GenomeThreader supports alignment of transcript evidence to multiple (nearly identical) loci. This is important for covering nearly identical tandemly-arrayed genes observed in P. patens[3, 10] and of other segmental duplications.

In V1.6 we have incorporated information on AS for the first time into the released genome annotation. The integration of AS into the annotation process leads to 3500 (10.8%) annotated loci with an average of 2.52 transcripts per locus and a maximum of 11 transcripts in V1.6. 1775 loci (51% of AS transcripts; 5.5% of all genes) have an altered coding sequence (CDS) due to alternative splicing resulting in 2380 distinct proteins. In contrast to A. thaliana, where the analysis of large-scale full-length cDNAs suggested most splicing events to occur outside of coding regions in the 5’-UTR [67], alternative splicing in moss seems to affect UTRs and CDS regions to a similar degree. 2948 alternative transcripts are due to alternative splicing in the UTRs of 1991 loci (56% of AS transcripts; 6.2% of all genes). Extensive alternative splicing of moss 5’-UTRs was observed previously in an individual study of the MDHAR genes [68].

Insights into the evolution of gene structures along the green lineage

Corresponding with its intermediary phylogenetic position, the median intron length of the average P. patens gene is intermediate (Table 3) between that of algae and flowering plants, implicating a trend towards shorter introns during the evolution of land plants. Whereas the median algal 5’-UTR length are in the range of Tracheophyta, the median moss 5’-UTR length is about 100 bp longer. Moss 5’-UTR regions more frequently comprise multiple exons than other plant genomes, with 40% of genes so annotated. Additionally, although there are substantially more predicted transcripts in O. sativa and more transcripts with 5’-UTR introns (8227), the total number of transcripts with multiple 5’-UTR exons in P. patens is nearly as high (7120) (p-value < 2E-238); Figure 2). This will likely increase as further high-resolution transcriptome data are obtained, since only 18,180 moss transcripts currently have annotated 5’UTRs compared with 31,793 transcripts in the rice annotation.

In A. thaliana introns in the 5’-UTRs tend to be longer than introns in the coding sequence and 3’-UTRs [70]. A plot of P. patens intron lengths (Figure 3) reveals a similar distribution. P. patens 5’-UTR introns are significantly longer than introns in the CDS (p-value = 0) and 3’-UTR (p-value <1.03E-193). In particular, the percentage of introns longer than 500 bp is much higher in the 5’-UTRs than in the CDS or in the 3’-UTRs. The mean 5’-UTR intron length for P. patens is 520 bp, whereas it is 264 bp for the CDS and 268 bp for 3’-UTR introns.

Furthermore, in A. thaliana the 5’-UTR introns are preferentially located close to the initiating ATG codon [70]. A comparative distribution plot of the proximate 5’-UTR intron distance to the start of transcription and translation reveals similar distributions in P. patens and A. thaliana (Figure 4). Certainly the closeness of 5’-UTR introns to the initiating ATG is more pronounced in A. thaliana (A. thaliana ~75% < 65 bp; P. patens ~50% < 65 bp; p-value < 9.7E-165). Chung et al. (2006) provided initial evidence for a positive effect of the presence of 5’-UTR introns on gene expression. These findings and previous reports suggest that properties of 5’-UTR introns influence gene expression level, regulation, translation and nonsense-mediated mRNA decay [71, 72]. As a consequence, the high percentage of genes with 5’-UTR introns, their unique length, and the fact that ~50% of P. patens loci with annotated alternative transcripts undergo alternative splicing in UTRs, suggests that the moss makes frequent use of this type of gene regulation.

Improved functional annotation of moss proteins

In summary, we improved the GO annotation of P. patens gene products in general and demonstrate that the annotation can serve as a solid basis for comparative and exploratory analyses. Compared to A. thaliana (TAIR10), where only 7% of all genes have no GO annotation, there are still more than 40% of all loci without any GO annotation for P. patens, and those with annotation are mainly inferred by electronic annotation (IEA). While a significant fraction of this 40% probably comprises orphan genes, special focus should be placed on the improvement of the functional annotation of the P. patens protein-coding genes to unravel its full potential as a reference model organism. To facilitate such efforts we developed the cosmoss.org community annotation interface to browse and alter the P. patens annotation (see section genonaut) discussed below.

Gene families

Gene families (clusters) were retrieved by clustering a multitude of sequences of Archaeplastida proteins (Additional file 5: Table A4) using OrthoMCL [79]. Parameters were optimized as described in Methods to target Archaeplastida gene families sensu strictu, i.e. families of genes that evolved by speciation and duplications after the divergence of the red/green lineage from a single gene in the last common ancestor. Multi-gene or superfamilies thus are split across multiple clusters. The number of P. patens loci in clusters is 32,733, while after the subtraction of the annotated number of 32,275 protein-coding loci, 458 fall into multiple clusters. Some of these cases are due to fragmentary or false structure predictions caused by fusion of two or more distinct gene loci into a single locus and will be resolved in future.

Strikingly, 48% of all P. patens loci are in Physcomitrella-only clusters (Figure 5), which independently supports an analysis based on protein domains where 52% of all P. patens genes have no Pfam domain [80]. Further filtering of these loci including BLASTP against GenPept and EST support indicates that ~22% (7169) of all loci in Physcomitrella-only clusters have no detectable homolog, while at least ~13% (4157) have no homolog but transcript evidence. These genes might represent true orphan genes, representing species- or lineage-specific adaptive innovations or non-protein-coding genes. The remainder falls into 8208 clusters with at least one other Viridiplantae protein (Figure 6). Thus, a prominent part of the P. patens gene complement could not be clustered with other plant proteins. In part, this could be explained by remaining fragmentary gene predictions and non-protein-coding genes. This is supported by the observation that the average transcript length of these P. patens-only singleton genes is less than 500 bp and the average number of exons is only two (total average is four exons, Table 3). However, at least 40% (4503) of these P. patens-only singleton genes have transcript evidence. Approximately two-thirds of all genes are in small clusters with less than three P. patens members. As already observed in the draft genome [10], several gene families are expanded by comparison to flowering plants. In total 7013 OrthoMCL clusters are shared by P. patens and A. thaliana containing 14,830 and 11,926 genes, respectively. Reported expansions in gene families [10], like LHCs (esp. chlorophyll a/b binding proteins), dyneins, histidine kinases and response regulators are also reflected by OrthoMCL cluster sizes (Figure 6). The majority of expanded families are in OrthoMCL clusters with more than five members in moss.

Analysis of intron-loss and gain in the green lineage

Similar to 5’-UTRs, where we observe a striking number of multi-exon regions in moss, the overall number of single-exon transcripts is remarkable. While Chlamydomonas and Volvox contain fewer than 10% single-exon transcripts (Table 3 and Additional file 4: Table A3), this fraction is on average ~19% in other land plants, while moss possesses 23.4%. This may be due either to fragmentary gene predictions and residual non-protein coding genes, or may reflect secondaryintron losses.It has been observed that introns and their positions are highly conserved during land plant evolution [81]. In A. thaliana and O. sativa intron losses outnumber intron gains [82]. Furthermore, there is evidence for secondary intron loss; e.g. a moss sedoheptulose-1,7-bisphosphatase (SBP) gene which lost six out of seven introns [81]. Such intron losses might account for the relatively high percentage of single-exon genes in P. patens. One suggested mechanism for intron loss involves the reverse transcription of an mRNA followed by the (partial) replacement of the genomic DNA copy by an intron-less cDNA via homologous recombination, called retrocopying [83]. An increased rate of intron loss might thus be facilitated by the extraordinarily high rate of DNA repair by homologous recombination in P. patens. This hypothesis is based on models proposing a prominent role of gene conversion and DNA repair in intron loss [84, 85].

About 3% (941) of the P. patens genes seem to have lost their introns entirely, i.e. are of putative retrocopy origin, which is in the range of the other land plants (Table 6). The extension of the analysis to Arabidopsis, Rice and Chlamydomonas supports the findings from the comparison of absolute numbers of single-exon genes between algae and land plants described in the previous sections. The alga has significantly less (0.7%) single-exon genes than the three land plants (3-4%) under study. One likely scenario is an increased activity of transposons resulting in a secondary, maybe more recent, intron gain in algae [69, 86]. This view is supported by the observation that intron positions are often not conserved between the two algae and the land plants (data not shown). Considering the comparable rate of intron loss in plant gene families, the fact that the total number of multi-exon transcripts in P. patens is similar to vascular plants, that only 60% of the 8979 single exon models are supported by expression evidence, and that more than half of these transcripts are shorter than 500 bp (which is less than half of the mean of P. patens transcripts, 1389 nt), leads us to conclude that a significant number of the predicted single exon genes represent fragmentary predictions, non-protein coding genes or pseudogenes.

Gene family size evolution in Viridiplantae

We compared gene family size distributions along the green lineage by looking at all clusters containing P. patens and at least one other Viridiplantae species (8208 clusters/families; Figure 6). This analysis revealed that protein families with two to six members are more frequent in P. patens than in the other Viridiplantae (Figure 7), supporting the previously hypothesized balanced retention of paralogs that act as pseudoalleles in a haploid-dominant land plant [3]. This is most pronounced in clusters with two members, which amount to 27% in P. patens while being substantially less frequent in the other green organisms. In total, 832 clusters point to an expansion in P. patens in comparison to the other Viridiplantae. 184 of these P. patens gene families are more than twice the size of their largest green lineage counterpart.

Among these expansions house-keeping and metabolic gene functions are most prominent, independently supporting analyses of paralogs retained after the proposed whole genome duplication event ~45 million years ago [87] and of the unique presence of identical tandemly-arrayed genes [10]. Prominent examples of this expanded category of genes are abundantly expressed components of multimeric protein complexes like the ribosome and proteasome. In total, 86 of the clusters represent the different abundantly expressed protein components of the ribosome. On average these harbor ~1.5 times more genes in P. patens than in A. thaliana. For the 36 clusters representing structural components of the proteasome we observed on average 1.6 times the gene complement of Arabidopsis in the moss. These expansions are also detectable in the functional comparison of the two species using GO enrichment analysis (Additional file 6: Table A5).

In addition, families encoding for smaller complexes and monomers are expanded, including the Light Harvesting Complex II (LHCII) major antenna, Ribulose-1, 5-bisphosphate carboxylase oxygenase (RuBisCO) small subunit, TOC12 (translocase of the outer chloroplast membrane 12) and components of the splicing and translation machineries.

The expansion of house-keeping and metabolic gene functions is also mirrored by the findings of the GO enrichment analysis comparing P. patens and A. thaliana (Table 5 and Additional file 6: Table A5), which revealed an increased complement of genes involved in translation, oxidation-reduction, electron transport, microtubule based movement, glycolysis, and ATP synthesis. Additionally, specific expansions occurred which might represent lineage-/species-specific adaptations (e.g. expansins, MIKC*-type MADS box transcription factors, late embryogenesis abundant (LEA) proteins, early response to dehydration (ERD) proteins, cationic peroxidases) and in enriched GO “biological process” annotations (e.g. chitin catabolic processes, the phosphoenolpyruvate-dependent sugar phosphotransferase system, cell wall macromolecule catabolic processes, phosphatidylinositol-mediated signaling, chromatin assembly or disassembly, cell redox homeostasis, ciliary or flagellar motility and DNA repair). Some of these enriched categories like “ciliary or flagellar motility”, which can be explained by the absence of flagellated sperm in flowering plants, confirm the findings of previous analyses [10, 69, 88]. The majority of the above listed categories and families (nine out of 13) represents true novel insights which will help to unravel the much-cited qualities of the mosses in coping with abiotic and biotic stressors like their unique ability to repair DNA damage [89–91].

Extension of the cosmoss.org resource to provide a permanent, central model organism database and annotation repository for P. patens

Protein-coding gene models were clustered to loci and the resulting locus definitions were used to derive information-rich locus identifiers (Cosmoss Gene ID; CGI; see Additional file 7: S2 and Additional file 8: Figure A2, and [92] for details).

With the development of genonaut [93], we have extended cosmoss.org by the capability to annotate Physcomitrella patens genes with regard to gene name, product name, description, and Gene Ontology (GO) terms. The interface which allows searching, browsing and editing of annotation is modular and can be extended to support additional (ontology) annotations as gene features. Traceability of annotations in terms of author and experimental evidence is crucial for quality assessment of information retrieval. Thus, the genonaut interface accepts the alteration of an existing gene description only if the source is specified. Integration of multiple sources into a unique abstraction layer is achieved by assigning unique and permanent Cosmoss Reference IDs (CRID) to author statements from all sources. The highest quality author statements are experimental evidence provided as references to peer-reviewed publications. The easiest way to achieve this is to provide a valid PubMed ID and the system automatically retrieves all relevant information from NCBI PubMed. If the source is a publication that is not tracked in PubMed, a custom reference can be created. If no publication is available, as the “weakest” possible evidence, a note in form of a text comment or web link describing the evidence is required. Besides the references the genonaut interface allows to link to other resources via database cross-references (Dbxref).

The gene products can be further annotated using GO terms. To allow convenient manual GOA, the genonaut interface assists the annotator by allowing the user to browse the appropriate GO namespace by keywords to assign the correct terms. In addition to the mandatory reference for a genonaut annotation, we have integrated the assignment of the GO evidence codes [75]. In this way the quality of each assigned GO term is directly discernible. As traceability is crucial for the maintenance of annotation quality, the genonaut system traces every annotation change using a history system. Thus it is possible to trace changes and possibly revert to a previous state if needed, but more importantly to comprehend the annotation history of every gene and annotation version.

Whereas the annotation browser capabilities are publicly available, the editor functions are restricted to registered cosmoss.org users (December 2012: 228 annotator accounts). The cosmoss.org curator team acts as a superior authority which supervises and validates the user provided annotations by direct personal communication.

Moreover, the genonaut interface provides a starting point to retrieve detailed annotation about P. patens protein-coding genes. Besides the possibility to search and edit the annotations and annotation history, the genonaut interface is linked to the sequence retrieval, genome browser and sequence viewer providing transcript and protein domain annotations.

To support the manual curation of gene structures we have integrated and adapted the Apollo structural gene annotation editor [94] for the cosmoss.org genome browser. Generated user_models (December 2012: 830 manually curated transcripts) are assigned CGIs (Additional files 7 and 8) from the user model namespace extended using the authors username to allow multiple versions per locus (e.g. Pp1s275_35U2 ➔ Pp1s275_35U2__zimmer.1).

P. patens reference genes

In total, 137 manually annotated and validated Physcomitrella patens gene structures are in the cosmoss.org genome browser track “Ppref genes”. Some of the genes are directly derived from published genes (GenBank) or provided by the scientific community, but the majority was extended or corrected manually using ESTs and FLcDNAs. In addition, in-house validated sequenced P. patens gene structures were added.

EuGène P. patens gene prediction process

EuGène [58] allows combining various type of evidence including e.g. ESTs, mate pair information, homologous sequences, existing gene predictions and splice site predictions. Our EuGène predictions for P. patens are based on a splice site prediction using SpliceMachine [57] trained on filtered P. patens EST alignments, coding, intronic and intergenic regions, homology evidence (A. thaliana homologs) and are filtered using transposon related sequences and gaps. The optimal parameters were determined on an independent P. patens reference genes set not used for training of EuGène. Two P. patens whole genome EuGène predictions went into generation of the release V1.6. The first contains 37,872 predicted loci and was restricted to generate UTR regions only if transcript evidence is available and the second contains UTR regions predicted ab initio (46,071 loci).

The first EuGène model predicts ~94% of all CDS exons in the reference genes correctly, whereas the second does so for ~95% (results are summarized in Table 3). Both predict 76% of all reference CDS without any error. That implies that the EuGène predictions perform well in predicting the exons but split several loci into two or more distinct genes. While manually inspecting these problematic loci we have noticed that gene models created by the JGI [10], which were not been selected for release V1.1, were often better than the selected model model, or could be used to overcome, or respectively complement, the EuGène predictions. As the method of choice to combine all available P. patens evidence and to further improve the protein-coding gene structures we have used EVidenceModeler [59].

EvidenceModeler (EVM) - weighted consensus gene model predictions

EVM (Haas et al. 2008) combines evidence from different sources into a consensus gene structure prediction. With the possibility to weight and unite the different evidence and optimize their combination, EVM utilizes the different sources by equating the drawbacks of individual sources but also boosting their strong points. For P. patens we have used EVM to find the optimal combination of PASA [60] transcript assemblies, EST alignments and five different whole genome protein-coding gene predictions. The process is also described in Additional file 3: Figure A1. The resulting models were subsequently subjected to PASA to model the UTR regions. As a consequence, all UTRs in release V1.6 are supported by transcript evidence. The utilization of EVM has enabled us to increase the prediction performance on the reference gene set (86.1% of all CDS and 97.3% of all CDS exons are correct; see Additional file 2: Table A2).

Additional gene structure predictions using EuGène

EuGène [58] version 3.4 was adapted and trained for P. patens on the basis of the Ppref genes set (mentioned above) and including a species-specific splice-site prediction and IMM (Interpolated Markov Model) models trained on intergenic, CDS and UTR regions. Gaps in the genomic sequence and repetitive regions, in particularly LTR retrotransposons [10], were masked for the training and predictions. As additional evidence, homologous protein sequences (Swissprot rel. 13.4 and Arabidopsis thaliana TAIR7 [34] homologs) and in particular EST alignments 360,974 from GenomeThreader [61], 118,243 from sim4 [95] and 97,373 from exonerate [96] were used. If available, we also provided EST mate pair information into EuGène (63,945 EST mate pairs). Two whole genome EuGène predictions were used for the consensus model approach leading to V1.6. The training input of these two models was the same except that one model was additionally trained with 5’- and 3’-UTR regions (EuGène MarkovIMM plugin).

Splice site prediction

P. patens EST alignments were loaded into a Bio::DB::Seqfeature database. The exclusion of alternative splicing and bad quality EST alignments led to a distinct, species-specific splice-site training set for SpliceMachine. Splice sites were only taken into account where GenomeThreader [61], sim4 [95] and exonerate [96] did exactly the same EST alignment. Two models, one donor (GT) and one acceptor (AG) model have been generated and have been used for EuGène training.

PASA – transcript assemblies

PASA [61] assemblies were performed as directed on the software homepage [97], with the following modifications: Per default PASA uses GMAP [62] for transcript alignments, however, our evaluation process (data not shown) reveals improved EST alignments using GenomeThreader. PASA offers the possibility to include alignments in GFF3 format via the --IMPORT_CUSTOM_ALIGNMENTS_GFF3 switch. The P. patens GenomeThreader EST alignments were converted into PASA compliant GFF3 format. PASA per default supports only one alignment per EST/transcript. Therefore with regard to duplicated genes (especially tandemly arrayed genes and (near-)identical genes) the corresponding EST alignments were renamed. E.g. the sequence ppsp14d22fl matches to 2 loci in the genome, so the corresponding sequence is duplicated and renamed: ppsp14d22fl and ppsp14d22fl_2. Transcript alignments with less than 90% EST length coverage or less than 95% alignment identity were discarded. The maximum allowed intron length was set to 20,000 nt, based on the longest observed intron supported by Sanger ESTs. Three cycles of PASA annotation loading, annotation comparison, and annotation updates were used to maximize the incorporation of transcript alignments into the transcript assemblies.

Consensus gene predictions – unfiltered V1.6 models

The P. patens JGI AllModels V1.1 (http://ftp.jgi-psf.org/pub/JGI_data/Physcomitrella_patens/v1.1/transcripts.Phypa1_1.AllModels.fasta.gz) and the cosmoss.org EuGène models, described in the previous section, and transcript alignments (PASA, exonerate and sim4) were combined and weighted with EVidenceModeler (EVM). EVM combines all evidence per locus into one consensus gene structure model. EVM was trained on the Ppref genes deducting reference genes used for the training of EuGène. Subsequently the UTR regions were modeled by P. patens PASA transcript alignment assemblies as described on the PASA homepage [97]. In this context models with alternative splicing evidence were generated and incorporated into V1.6.

Functional gene annotation

The unfiltered V1.6 gene models and the corresponding predicted protein sequences, respectively, were subjected to BLAST2GO [73]. The initial BLASTP (e-value cut-off: 1E-4) search was performed against a P. patens V1.1 subtracted GenPept release 172.0 to allow an unbiased functional annotation, independent from the previous entirely automatic V1.1 GOA. For the BLAST2GO annotation step the minimum coverage between a hit and its HSP was set to 40%. The validation step as well as the integration of InterProScan V4.5 (InterPro release v22.0) based GO annotation was used to generate the GOA for the P. patens V1.6 proteins.

GO term annotation was extended by experimental evidence, various subcellular target predictions and homology-based methods using the pred2GOA method described below. Existing functional annotations like gene names, description lines, GO terms and KEGG EC numbers and KO terms were collected from JGI and Kyoto Encyclopedia of Genes and Genomes [98] and combined into a non-redundant database. Existing experimental evidence for subcellular localization was also manually integrated as GOA if available. Additional functional annotations were created using the homology-based methods BLAST2GO, IPRScan [74] and KAAS [99]. These two steps resulted in GO terms with IEA evidence codes. Subcellular localization of the protein sequences was predicted considering the individual gene’s full-length status using several tools: TAPScan [100], MultiLoc [101], WolfPSORT [102], TargetP [103], ChloroP [104], SignalP [105], Prosite KDEL [106], HMMTOP [107] and MEMSAT3 [108]. To combine these predictions with the existing GOA, we developed the algorithm pred2GOA which allows weighted integration of GOA from multiple sources. Resulting predictions were translated into cellular component GO terms and compared to the existing GOA at the GO slim (plants) level. The assignment of GO terms is based on a weighted majority rule consensus. If the underlying gene did not have an annotated UTR and a start codon, at least one of the predictors’ predictions had to be based on more than the N-terminal region of the protein. Resulting GO term assignments were reviewed for consistency and defined as “Inferred from Sequence or Structural Similarity” (ISS). The resulting GO annotation was mapped to GO slim terms using the Blast2GO internal mapping function using the “goslim_plant.obo” ontology subset.

GO enrichment analysis

The enrichment analyses were performed using the Bioconductor package topGO [109]. We used and compared results from both the classical Fisher’s Exact Test and the topGO algorithm “weight01” for test statistics with a p-value cut-off 0.05. The Arabidopsis GOA was downloaded from the TAIR ftp server (TAIR10; http://ftp.arabidopsis.org/Ontologies/Gene_Ontology/ATH_GO_GOSLIM.txt).

Protein clusters and representative gene model selection

Although certainly useful in single gene analysis and the study of AS, the use of all splice-variants in large-scale comparative analysis introduces undesirable complications. Therefore gene catalogues are usually reduced to one representative isoform per locus prior to large-scale analyses. By convention, the representative model is the variant with the lowest splice variant index. A representative gene model per locus was selected by the following procedure: At the first step, all P. patens V1.6 proteins together with the proteins of various sequenced Archeaplastida (Additional file 5: Table A4) were clustered using OrthoMCL V.2 [79]. The inflation parameter was optimized using a set of reference gene families (e.g. [10, 15]). The inflation value was set to 1.3. The resultant clusters were aligned with MAFFT ginsi v6.712b [110] and subsequently the uncorrected pairwise distances were calculated with distmat (EMBOSS 6.1.0; [111]). The distmat matrices were used to select the Physcomitrella model as representative for a locus that has the least of all substitutions (minimum distmat distance) to any non-Physcomitrella cluster member. Clusters with only Physcomitrella sequences were subjected to a BLASTP search against GenPept (release 172.0). In this case the representative per locus was set to the Physcomitrella sequence which covers its closest GenPept hit, in terms of alignment length, best. If still no clear representative could be determined, the model with the longest transcript length was chosen as the representative. Due to the fact that the clustering is based on sequence information alone, not all splice variants per locus were grouped into the same cluster in some cases. Thus, the number of loci in clusters slightly exceeds the number of physical loci.

RNA genes

Pre-tRNA genes were predicted by combining results from tRNAscan-SE 1.21 [112] and ARAGORN 1.1 [113]. rRNA loci were predicted using RNAmmer-1.2 [114] and extended by mapping the available SILVA Physcomitrella rRNAs [115] to the genome with BLASTN.

Non-protein coding loci and gene families where determined using Infernal [116] with the RFAM (version 8.1) [117] covariance models and by integration of the miRBase 18 [51]Physcomitrella miRNA classifications and annotations.

P. patens V1.6 protein-coding genes filtering

On the basis of a Bio::DB::SeqFeature database [118] the V1.6 gene models were filtered against the annotated LTR-retrotransposons [10], Repbase (RELEASE 20080801; [119]) using RepeatMasker v3.26 [120] and the non-coding RNA described in the previous section.

Codon usage

The effective number of codons was calculated with CodonW 1.4.4 [121].

Statistics

Fisher’s exact tests and Wilcoxon tests were performed with R 3.0.0. The p-values were corrected for multiple testing using fdr.

Availability of supporting data

The complete annotation and sequence information comprising: Gene structure releases in GFF3 format: gene, transcript, CDS, protein, UTR sequences in FASTA format; Mappings/lookup tables, GO annotations in GAF format are accessible via the cosmoss.org download section [122].