Comparative genomics of Physcomitrella patens gametophytic transcriptome and Arabidopsis thaliana: Implication for land plant evolution

Plant Materials, RNA Preparation, and Construction of Full-Length cDNA Libraries.P. patens (Hedw.) Bruch and Schimp subspecies patens Tan (20), collected in Gransden Wood, Huntingdonshire, U.K. (21), was used as the wild-type strain. The protonemata were ground with the Polytron (Kinematica, Littau, Switzerland), and inoculated into BCDATG medium (22), BCD medium (22) that contained 1.0 mM CaCl₂ and 1.0 μM naphthalene acetic acid (NAA; Sigma), or BCD medium that contained 1.0 mM CaCl₂ and 0.50 μM 6-benzylaminopurine (BA; Sigma) for the untreated, NAA-treated, and BA-treated specimens, respectively, at 25°C under continuous light, and the tissues were harvested at 13–14, 8–11, and 8–13 days, respectively. At harvest, the untreated tissues contained chloronemata and young gametophores with two to five leaves, the NAA-treated tissues contained chloronemata, caulonemata, and rhizoid-like protonemata, and the BA-treated tissues contained chloronemata, caulonemata, and malformed buds.

Total RNA was extracted from each tissue sample as described (23), and poly(A)⁺ RNA was purified with oligo(dT) magnetic beads (DynaBeads; Dynal A.S., Oslo). The cDNA was synthesized by using trehalose-thermoactivated reverse transcriptase (24), and full-length cDNA was recovered by using the biotinylated CAP trapper method (25, 26). The λ-full-length cDNA vector (27) and the single-strand linker ligation method (28) were used in the construction of the cDNA libraries. One round of normalization (29) was performed to construct the NAA- and BA-treated cDNA libraries to reduce the frequency of highly expressed mRNAs in the libraries.

DNA Sequencing. In the preliminary analyses, 908 clones were chosen from the untreated library, and the DNA was sequenced from the 5′-end. For mass sequence analyses, the clones from the untreated, NAA-treated, and BA-treated libraries were arrayed in 35, 50, and 48 384-well-plates, respectively, and the clones were sequenced at both ends. Low-quality sequences were removed based on PHRED quality score (30) so that average score is at least 20. The sequence data from six plates were inconsistent in terms of their 5′ and 3′ end sequences because of mishandling; therefore, the reverse sequences were renamed and treated as independent sequences during later analyses.

Sequence Data Sets. In addition to the 85,191 ESTs sequenced in this study, all of the P. patens mRNA entries in GenBank release 131.0 (excluding our ESTs) were extracted. ESTs with the annotation of “High quality sequence stop” were trimmed to the length specified by the annotation. All of the A. thaliana protein sequences were retrieved from the RefSeq database (31) and used in the A. thaliana protein data set. The Synechocystis sp. PCC 6803, Saccharomyces cerevisiae, and Drosophila melanogaster coding sequences and translated sequence data sets were obtained from the following locations.

The nonredundant protein sequence data set was obtained from ftp://ftp.ncbi.nih.gov/blast/db/nr.Z.

Identification of Taxonomic Positions. The taxonomic classification of the organisms followed that of the National Center for Biotechnology Information taxonomy database (32). The Gen-info Identifier (GI) number for the taxonomy ID (tax ID) table was obtained from ftp://ftp.ncbi.nih.gov/pub/taxonomy/gi_taxid_prot.dmp.gz, and the tree structure was according to “nodes.dmp” in ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz. The tree was traversed from a node that was identified by the tax ID to the root until a specified node was reached. When the root was reached, the taxonomic position was unidentified and thus investigated manually.

Statistical Test of Transcript Representation in Different Libraries. The number of clones in each library was counted for each putative transcript. We used the null hypothesis that a putative transcript is equally represented in the three libraries and therefore the proportion of clones belonging to the putative transcript is constant in the three libraries. We used the G test of independence (33).

Library Quality. Full-length cDNA clones are important, not only for gene annotation and the determination of transcriptional start sites, but also for functional analyses (34, 35). Recently, methods have been developed that allow preferential cloning of cDNA that corresponds to full-length mRNAs with 5′-end-proximal cap structures (reviewed in ref. 36). These methods have been used in large-scale analyses of transcripts from human (37), mouse (38, 39), fruit fly (40), rice (39), and A. thaliana (41, 42).

We constructed three full-length enriched cDNA libraries from the gametophytic tissues of P. patens that were NAA-treated, BA-treated, or untreated, and we selected 19,200, 18,432, and 14,348 clones, respectively, from the libraries. The clones were sequenced from both ends. In total, 4.78 × 10⁷ nucleotides of 85,191 sequences remained after the elimination of low-quality sequences and contaminated clones showing nearly complete match to Escherichia coli or Dictyostelium discoideum sequence.

We analyzed the proportion of full-length cDNA clones that contained mRNA start sites. Because transcription usually starts at a purine base (43), we examined the base composition of the first and subsequent nucleotides at the 5′ end of each 5′ sequence. The first nucleotide showed a strong bias for purine (81% G or A; Table 1), whereas 42–44% of the 2nd to 20th nucleotides were purines. If the cDNAs were not full-length and did not contain the first nucleotide of the mRNA, then the base composition of the first nucleotide would have been identical to that of nucleotides in subsequent positions. Although terminal transferase activity of the reverse transcriptase, adding C to the 3′ end of cDNA (44), may have contributed partly to over representation of G, this does not account for the increase of A. Therefore, the deviation of first nucleotide base composition from that of nucleotides at subsequent positions suggests that a high proportion of the cDNA clones are full-length. In our 5′ end sequence data (Table 1), the T/C ratio at the first position (0.65) differed from that at the second (2.2) and subsequent (1.1–2.6) positions. This finding implies that a significant proportion of the 5′-end-proximal pyrimidine bases also represent mRNA start sites.

To estimate the proportion of clones that contained complete coding sequences, we examined two P. patens tubulin-α genes (Y. Sato, T.F., and M.H., unpublished data) and all of the 38 P. patens mRNA deposited in GenBank (release 130.0) that were annotated as containing complete coding sequences. ESTs corresponding to the genes were identified with blastn (45) searches and then tested whether they had 5′ or 3′ UTR. The search revealed that 31 genes had corresponding clones in our EST library (Table 4, which is published as supporting information on the PNAS web site, www.pnas.org). In the 5′ end sequences, 95.4% (456 of 478) of the clones contained 5′ UTR, and 97.8% (480 of 491) of the 3′ end sequences contained 3′ UTR (Table 4). If we assume that 5′ and 3′ end truncation occurred independently, the proportion of clones that had both start and stop codons can be calculated as the product of the proportion of clones complete at the 5′ and 3′ ends. Therefore, we estimated that 93% (95.4% × 97.8%) of the clones contained complete coding sequences.

Assembling the EST Data into Contigs and Estimating the Gene Number.P. patens mRNA sequence data from other sources, such as the ESTs from the Physcomitrella EST Programme, were obtained from GenBank (release 131.0) and pooled with our data. After removing low-quality sequences, which were identified by their annotations, and vector and linker sequences, which were excluded based on sequence similarity searches, a data set of 5.4 × 10⁷ total nucleotides in 102,553 sequences was obtained. The data were assembled with the CAP3 Sequence Assembly program (46). We obtained 22,885 contigs in total. Of these, 20,589 were represented in our library, whereas the remaining 2,296 contigs arose from other sources. Because our ESTs were obtained by sequencing each clone from the 5′- and 3′-ends, we can assign 5′- and 3′-contigs to the original clone. In total, 814 of the 20,589 contigs had >5% ambiguity (N), and were represented by only one sequence (singlet). These contigs were not used for pairing, because they may not represent independent transcripts. By using the data on the 5′- and 3′-ends, the remaining 19,775 contigs were organized into 13,593 putative transcripts. We assigned a unique identifier (of the form Pnnnnnn, where n is a digit) to each of the 13,593 putative transcripts. Each of the 2,296 contigs from other sources, with the exception of six contaminating sequences from bacteriophage λ, was regarded as an independent putative transcript. The addition of these 2,290 putative transcripts to the 13,593 putative transcripts reveals that at least 15,883 putative transcripts are expressed in P. patens.

Although >100,000 sequences were determined for the 22,885 contigs, 9,417 of the 22,885 contigs were composed of only one sequence, i.e., either an EST or a deposited mRNA sequence. This finding indicates that the EST library is not saturated, and the sequencing of more clones will probably reveal new genes. In addition, gametophytes with gametangia and sporophytes were not sampled. Therefore, the total number of P. patens transcripts probably exceeds 15,883.

Some Genes Are Represented Differently in the Three Libraries. Because we sequenced three libraries, we were able to identify genes that were represented differently in each library. We used the G test (33) to assess whether putative transcripts were over- or underrepresented in the auxin- and cytokinin-treated libraries. Putative transcripts of this type contained genes that were regulated by auxin and cytokinin, as well as genes that were related to differentiation and development induced by the auxin and cytokinin treatments.

Seven PpIAA1 cDNA (14) clones were found in the auxin-treated library and one clone was found in the cytokinin-treated library. In this case, the null hypothesis of equal representation in the three libraries was rejected at the 1% level. Induction of the PpIAA1 gene by auxin is consistent with the reported RNA gel blot analysis (14).

Similarly, all of the putative transcripts were tested for over- or under-representation in each library. Sixty and sixty-eight putative transcripts were over-represented at the 0.1% level in the auxin- and cytokinin-treated libraries, respectively (Tables 5 and 6, which are published as supporting information on the PNAS web site). One and twenty-one putative transcripts were under-represented in the auxin- and cytokinin-treated libraries, respectively (Table 7, which is published as supporting information on the PNAS web site).

Of the putative transcripts that were overrepresented in the cytokinin-treated library, two (P013130 and P006794) showed high-level similarity to cytokinin oxidase, which functions in cytokinin degradation (47). The promoter of the cytokinin oxidase gene DSCKX1 of Dendrobium cv. Sonia was shown to be cytokinin-inducible (48). Overrepresentation of the cytokinin oxidase gene in the cytokinin-treated library is consistent with induction by cytokinin, and suggests that mosses and vascular plants share common mechanisms for the regulation of cytokinin levels. The receptor-like kinases (P007265, P008576, P001848, and P006061) and the ethylene responsive element binding protein (EREBP)-like transcription factors (P012262 and P001581) were also overrepresented in the cytokinin-treated library, though the functions of these factors are unknown.

Many of the A. thaliana Genes Have Homologues in the Moss Gametophytic Transcriptome. To analyze the proportion of A. thaliana genes that were expressed in the moss gametophytes, all 26,178 amino acid sequences encoded by A. thaliana genes were obtained from the RefSeq database (31) at the National Center for Biotechnology Information and used as the query sequences in a similarity search with tblastn (45) to the moss data set. The peak occurred at 0.1 < E ≤ 10 and the tail appeared at E ≤ 1e-3 (Fig. 2, which is published as supporting information on the PNAS web site). If we regard the P. patens contigs with E ≤ 1e-3 as homologues of A. thaliana query sequences, 17,382 of 26,178 (66.4%) of the A. thaliana genes have homologues in P. patens. Random matches with an E value ≤ 1e-3 would be expected to occur only 26.2 (26,178 × 0.001) times. Matches occurred at a much higher frequency than this, and most were considered to be caused by homology.

tblastn similarity searches were also performed by using the amino acid sequences of the 26,178 A. thaliana genes as query sequences and the sequence data sets from budding yeast (S. cerevisiae), fruit fly (D. melanogaster), and cyanobacterium (Synechocystis sp. PCC 6803); putative homologues (E ≤ 1e-3) were discovered at frequencies of 40.5%, 47.1%, and 21.3%, respectively. The overlaps of A. thaliana genes with homologues in different organisms are summarized in the Venn diagram in Fig. 1. Of the A. thaliana genes, 5,039 showed significant similarities to P. patens genes, but not to genes of the other three organisms.

Sixteen classes of plant-specific transcription factors have been reported in the A. thaliana genome (10). Of these, the MYB R2R3 (49), homeodomain-leucine zipper (50), and ARF-AUX/ IAA (14) genes have been reported, and these genes were also found in the moss EST data. The NAC-like, zinc finger-DOF, zinc finger-WRKY, AP2/EREBP, ABI3/VP1-like and GRAS-like transcription factors were found in the EST data set. The putative transcript P011588 showed high-level similarity to Mybst1, which is a Myb single-repeat protein, which suggests that P011588 also encodes a Myb single-repeat protein, though the complete sequence is required to verify that P011588 encodes only a single Myb domain. The remaining six classes of transcription factors (cycloidea-like, YABBY, EIN3, RAV-like, squamosa-binding protein, and GT1-like) were not identified in the P. patens EST data set.

Some gene families, such as the REL-homology region protein, nuclear steroid receptors, forkhead-winged helix, and POU were found in metazoans, but were absent from the A. thaliana genome (10) and P. patens EST data set.

Because 17,382 (66.4%) of the predicted A. thaliana genes had related sequences in the P. patens ESTs and the ESTs are not saturated, the number of genes with P. patens homologues should be higher. This observation indicates that the P. patens gametophytic transcriptome is substantially similar to the A. thaliana genome. Although the A. thaliana genome sequence has been largely determined, the functions of many of the predicted genes are still unknown. A number of A. thaliana genes have homologues in P. patens, but not in the yeast or fruit fly. Because gene targeting is feasible in P. patens, functional analyses of the homologues in P. patens should increase our understanding of gene function in flowering plants.

P. patens Expresses Genes That Are Not Present in the A. thaliana Genome. We assessed whether the closest relatives of moss genes existed in vascular plants or in other organisms. The putative moss transcripts were used as the query sequence in a blastx (45) search against the A. thaliana protein data set and the nonredundant amino acid sequence (nr) data set in GenBank. In total, 9,942 (62.6%) of the 15,883 putative transcripts had significant similarities (E ≤ 1e-3) to A. thaliana proteins, and 10,757 (67.7%) had significant similarities to the nr data set. This is substantially higher than the estimation by Rensing et al. (19) that only half of the P. patens genes have known homologue in plants, fungi, animals, or bacteria, despite the fact that we used only gametophytic transcripts. This discrepancy is perhaps because we connected 5′ and 3′ contigs into a single putative transcript, which was used for the search. When we searched with the contigs as query sequences, half of the contigs showed no significant similarity to known proteins, which is consistent to the result by Rensing et al. (19).

We sorted the 10,757 putative transcripts into groups based on the organisms that had the highest levels of similarity to the moss transcripts. Of these, 9,907 transcripts had highest levels of similarity to genes in vascular plants, whereas 850 showed highest levels of similarity to genes in other organisms (Table 2 and Tables 8 and 9, which are published as supporting information on the PNAS web site). Some of the 850 putative transcripts showed similarities to both vascular plants and other organisms, with low levels of difference. Because the significance of these differences could not be ascertained, we selected 300 of 850 putative transcripts with E values that were at least 1e7-fold lower than the E values for best hits with vascular plants. These are candidates for transcripts that are absent in vascular plants but present in mosses. We termed these 300 putative transcripts MSTs (moss transcripts absent in vascular plants). Thirteen of the 300 MSTs were found in liverworts, which are also bryophytes.

The 300 MSTs probably correspond to new genes in land plants. Genome sequencing projects with several vascular plants are underway. Therefore, genes with high levels of similarity to some of the 300 MSTs may be found in future genomic analyses. Nonetheless, we believe that the 300 MSTs may be attributed primarily to gene loss or rapid evolution in vascular plant lineages, or in some cases, to horizontal gene transfer to P. patens from other organisms. Sequence conservation among the mosses and other organisms implies that the genes are important in mosses, whereas genes with lower or null levels of similarity in the A. thaliana genome either have changed function or are not required in flowering plants, which likely include genes that are involved in moss-specific processes.

Characterization of Moss Transcripts That Are Not Present in Vascular Plants. The 9,907 putative transcripts with the highest levels of similarity to vascular plants corresponded to 6,773 genes in the nonredundant data set, and the 300 MSTs corresponded to 244 genes. The gene ontology (GO) (51) categories of biological processes were assigned to 2,946 of the 6,773 genes and 116 of the 244 genes, respectively, based on INTERPRO (52) entries that matched the amino acid sequences of the 6,773 and 244 genes. The remaining 3,827 and 128 genes did not match INTERPRO entries. The frequencies of the GO categories at the fourth level are summarized in Table 3. The composition of the annotations was generally similar between the 2,946 and 116 genes, in that the “metabolism” category appeared with the highest frequency; the percentage of genes in this category was lower in the 116-gene group (68.0%) than in the 2,946-gene group (77.7%). The category “metabolism” was divided into subcategories; of these the “protein metabolism” and “biosynthesis” subcategories had fewer representatives among the 116-gene group (12% and 7% of metabolism, respectively) than among the 2,946-gene group (20% and 14%, respectively; Table 10, which is published as supporting information on the PNAS web site). These results suggest that genes that are involved in protein metabolism and biosynthesis are well conserved between mosses and vascular plants. By contrast, the group of 116 genes had a higher proportion of genes that were associated with transport and signal transduction than the group of 2,946 genes (Table 3). This finding may be related to the fact that most cells of mosses, especially the protonemata, are in direct contact with the outside environment.

Examples of Moss Transcripts That Are Absent in Vascular Plants.P. patens is known for its extraordinarily high ratio of homologous recombination to nonhomologous recombination. Some of the genes that had the highest levels of similarity to genes from organisms other than vascular plants were related to DNA damage repair, i.e., the P011933 gene for the DNA repair hydrolase, the P002718 gene for RecA, and the P009868 gene for Eso1. P009868 showed high-level similarity to the Eso1 gene of the fission yeast Schizosaccharomyces pombe, with a score of 120 bits and an E value of 7e-27. Eso1 is essential for sister chromatid cohesion (53), which is important for the efficient homologous recombinational repair of DNA double-strand breaks, because the intact sister chromatid paired to the damaged chromatid serves as a repair template (54, 55). Eso1 is a fusion of DNA polymerase η, which functions in DNA damage repair, and the sister chromatid cohesion protein (53). Although a high level of similarity was found within the DNA polymerase η domain, it remains to be seen whether the moss gene has the sister chromatid cohesion domain, because we obtained only the end sequences. DNA damage repair may be more important during the haploid stage than during the diploid stage, because only one copy of each gene exists in the haploid, whereas an extra copy exists in the diploid cell. Analysis of these genes may unveil the reason why P. patens has such a high rate of homologous recombination.

The putative transcripts P001256 and P005876 were highly similar to the chlorophyll a/b-binding protein (cab) homolog LI818r-1 of Chlamydomonas reinhardtii. The cab homologue LI818 was first identified in Chlamydomonas eugametos as a light-inducible mRNA that encoded a new class of the cab gene family (56). This class of cab homologues has not been reported in other land plants, and may have been lost during the evolution of the vascular plant lineage.

P006275 is an example of a gene that was probably transferred from bacteria, because it showed highest similarity to the glutatione S-transferase of Mesorhizobium loti (NP_105706), with a score of 251 bits and an E value of 9e-66. P006275 also had similarities with a number of glutathione S-transferases from other bacteria. The eukaryotic gene with highest similarity to P006275 was the glutatione S-transferase gene of Emericella nidulans, with 158 bits and 1e-37, whereas the closest gene in green plants was the glutatione S-transferase gene of Glycine max, with 29 bits and 1e-10. Clones of P006275 were found in all three of our libraries, and also in two libraries from the Physcomitrella EST Programme project, which indicates that this sequence is not a contaminant. The absence of a closely related homologue in plants was further tested in a tblastn search that matched the amino acid sequence of P006275 with the EST data set. The best non-P. patens hit was an EST from Porphyra yezoensis, which is a red alga, with 108 bits and an E value of 2e-22; no ESTs from green plant species were found, despite the existence of extensive EST libraries from several seed plants. Phylogenetic analysis (Fig. 3, which is published as supporting information on the PNAS web site) also supports a bacterial origin for P006275.

Genes of Sporophyte- and Gametophyte-Dominant Plants. Based on the phylogenetic analyses, the closest relatives of the land plants are the freshwater green algae of Charophyceae (57). The single diploid cells of charophycean algae are zygotes, which undergo meiosis to produce four haploid gametophytic cells. Within the land plant lineage, bryophytes diverged at the most basal position. Both charophycean algae and bryophytes have gametophyte-dominant life cycles, which suggests that the ancestors of land plants, from which sporophyte-dominant vascular plants evolved, had gametophyte-dominant life cycles. Because the multicellular sporophytes probably evolved during land plant history from the gametophyte-dominant plants in the absence of multicellular diploid plant bodies, it seems reasonable to speculate that the genes that function in sporophytes were recruited from genes that functioned in the gametophytes of the common ancestor of mosses and flowering plants. However, there is no direct evidence for this evolutionary scenario, because gene expression in gametophytes of land plants is not well characterized, even in model plants such as A. thaliana, in which almost the whole genome has been sequenced. The general similarity between P. patens gametophyte ESTs and the A. thaliana genome supports the evolutionary scheme that proposes the recruitment of gametophytic genes into the sporophyte generation during land plant history, though we do not have an accurate estimation of the number of A. thaliana genes that are expressed in sporophyte generation. Analysis of the genes that are expressed in A. thaliana gametophytes will provide further insights into the evolution of the gametophyte and sporophyte generations.

The gametophyte transcriptome of P. patens is similar to the A. thaliana genome. Our findings that >66% of A. thaliana genes have homologues in P. patens gametophytes, and that >90% of the most closely related homologues of P. patens gametophytic transcripts occur in vascular plants, suggest that gemetophytes and sporophytes use similar gene sets. On the other hand, 850 putative transcripts of P. patens showed lower levels of similarity to genes of A. thaliana and vascular plants than to genes in other organisms, thereby indicating that these genes are probably associated with moss-specific features, such as morphology, life cycle, responses to the environment, metabolism, and other aspects. Analyses of these genes will give further evolutionary insights and new genetic resources for agricultural purposes.

The functions of many A. thaliana genes have not yet been revealed. Because gene targeting is feasible in P. patens, analysis of related genes in P. patens will be a good way of achieving an understanding of the function of plant-specific genes. Our full-length cDNA collection is a valuable resource for functional genomics in plants and is distributed by the RIKEN Bio Resource Center (www.brc.riken.go.jp/lab/epd/Eng. The data can be accessed through internet at http://moss.nibb.ac.jp/.

Abbreviations: MST, moss transcripts absent in vascular plants; NAA, naphthalene acetic acid; BA, 6-benzylaminopurine.

We thank K. Oishi, S. Haga, T. Morishita, S. Miura, F. Ohta, H. Hayashi, S. Nishizaka, H. Nomoto, and M. Sano for EST sequencing. We thank Drs. A. Kitayama and N. Ueno for providing equipment and valuable discussions. We also thank M. Naruse, Y. Bitoh, and Y. Dodo for technical assistance. This work was mainly supported by a Grant-in-Aid for Scientific Research on Priority Areas (C) “Genome Science” from the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT) (to M.H., Y.K., M.S., and K.S.), Grant-in Aid for Scientific Research on Priority Areas (Grants 14036229 and 13044001) from MEXT (to M.H. and T.F.), a Grant-in-Aid from MEXT (to M.H. and T.F.), and Japan Society for the Promotion of Science (JSPS) Grant RFTF00L0162 (to M.H. and T.F.). Construction of the full-length library was also supported by a grant for Genome Research from RIKEN, the Program for Promotion of Basic Research Activities for Innovative Biosciences, the Special Coordination Fund of the Science and Technology Agency, and a Grant-in-Aid from MEXT (to K.S.). T.N. is a JSPS Fellow.