Source: https://www.nature.com/articles/nature16548?error=cookies_not_supported&code=183bdadd-f2e3-4bde-9489-63bcfe317d42
Timestamp: 2019-04-19 08:53:12+00:00

Document:
Seagrasses colonized the sea1 on at least three independent occasions to form the basis of one of the most productive and widespread coastal ecosystems on the planet2. Here we report the genome of Zostera marina (L.), the first, to our knowledge, marine angiosperm to be fully sequenced. This reveals unique insights into the genomic losses and gains involved in achieving the structural and physiological adaptations required for its marine lifestyle, arguably the most severe habitat shift ever accomplished by flowering plants. Key angiosperm innovations that were lost include the entire repertoire of stomatal genes3, genes involved in the synthesis of terpenoids and ethylene signalling, and genes for ultraviolet protection and phytochromes for far-red sensing. Seagrasses have also regained functions enabling them to adjust to full salinity. Their cell walls contain all of the polysaccharides typical of land plants, but also contain polyanionic, low-methylated pectins and sulfated galactans, a feature shared with the cell walls of all macroalgae4 and that is important for ion homoeostasis, nutrient uptake and O2/CO2 exchange through leaf epidermal cells. The Z. marina genome resource will markedly advance a wide range of functional ecological studies from adaptation of marine ecosystems under climate warming5,6, to unravelling the mechanisms of osmoregulation under high salinities that may further inform our understanding of the evolution of salt tolerance in crop plants7.
Seagrasses are a polyphyletic assemblage of basal monocots belonging to four families in the Alismatales1,2 (Supplementary Note 1.1 and Supplementary Fig. 1.1). As a functional group, they provide the foundation of highly productive ecosystems present along the coasts of all continents except Antarctica, where they rival tropical rain forests and coral reefs in ecosystem services8,9. In colonizing sedimentary shorelines of the world’s ocean, seagrasses found a vast new habitat free of terrestrial competitors and insect pests but had to adapt to cope with new structural and physiological challenges related to full marine conditions.
Zostera marina (Zosteraceae), or eelgrass (Fig. 1), is the most widespread species throughout the temperate northern hemisphere of the Pacific and Atlantic10. A clone of Z. marina was sequenced from the Archipelago Sea, southwest Finland, using a combination of fosmid-ends and whole-genome shotgun (WGS) approaches (Methods, Supplementary Note 2). The 202.3 Mb Z. marina genome encodes 20,450 protein-coding genes, 86.6% of which (17,511 genes, Supplementary Note 3.1) are supported by transcriptome data from leaves, roots and flowers (Extended Data Fig. 1, Supplementary Notes 3.2–3.3 and Supplementary Data 1–3). Genes are located in numerous gene-dense islands separated by stretches of repeat elements accounting for 63% of the non-gapped assembly (Extended Data Fig. 2, Supplementary Note 3.1) as compared to only 13% in the only other sequenced alismatid, the freshwater duckweek, Spirodela polyrhiza (Alismatales, Araceae)11. Gypsy-type (32%) and Copia-type (20%) transposable elements contribute to most of the repetitive DNA. Sequence divergence analysis suggests that the genome retains copies from two distinct periods of invasion by Copia elements, but only one period for Gypsy elements (Extended Data Fig. 3a–c). Genes gained by Z. marina (‘accessory’) are located closer to transposable elements than to conserved (‘single copy’) genes (Fisher’s exact test, P < 0.0001) indicating that transposable elements may have played a role in genic adaptation.
Figure 1: Zostera marina and phylogenetic tree showing gene family expansion/contraction analysis compared with 13 representatives of the Viridiplantae.
a, Gains and losses are indicated along branches and nodes. The number of gene families, orphans (single-copy gene families) and number of predicted genes is indicated next to each species. Background colours (top to bottom) are Alismatales, other monocots, dicots, mosses/algae b, Typical Zostera marina meadow, Archipelago Sea, southwest Finland (photo by C.B.).
We identified 36 conserved microRNAs with high confidence and their predicted targets (Supplementary Note 3.4, Supplementary Data 4 and 5). A novel variant of miR528 (not present in Spirodela) was found to be the only member of this miRNA family, and demonstrates that this conserved miRNA is the only one ancestral to the entire monocot lineage. Most likely, Z. marina did not take part in the subsequent birth of miRNAs that are common to several other monocots12; nor did it experience or retain traces of prominent miRNA duplications.
Analysis of synonymous substitutions per synonymous site (KS) age distributions indicates that Z. marina carries the remnants of an independent, ancient whole-genome duplication (WGD) event (Fig. 2a, Supplementary Note 4.1)13. Duplicated segments account for ~9% of the Z. marina genome, probably an underestimate due to the fragmented nature of the assembly. Zostera and Spirodela diverged somewhere between 135 and 107 million years ago (Mya)14 and phylogenomic dating13 of the Z. marina WGD suggests that it occurred 72–64 Mya (Fig. 2b), thus independently from the two WGDs reported for S. polyrhiza11. This timeframe coincides with the initial diversification of a freshwater clade that includes three of the four families of seagrasses (Supplementary Table 1.1) and with the Cretaceous–Palaeogene (K–Pg) extinction event (Fig. 2c), which provided new ecological opportunities and may have triggered seagrass adaptive radiations.
Figure 2: Ancient whole-genome duplication (WGD).
a, KS-based age distribution of the whole Z. marina paranome. The x axis shows the synonymous distance until a KS cut-off of 2, in bins of 0.04, containing the KS values that were used for mixture modelling (excluding those with a KS ≤ 0.1). The component of the Gaussian mixture model plotted in red (as identified by EMMIX) corresponds to a WGD feature based on the SiZer analysis (other components are shown in black). The transition from the blue to the red at a KS of ~0.8 in the SiZer panel (below) indicates a change in the distribution and therefore provides evidence for an ancient WGD (Supplementary Table 4.1, Supplementary Fig. 4.1). b, Absolute age distribution obtained by phylogenomic dating of Z. marina paralogues. The solid black line represents the kernel density estimate (KDE) of the dated paralogues and the vertical dashed black line represents its peak, used as the consensus WGD age estimate, at 67 Mya. Grey lines represent the density estimates from 2,500 bootstrap replicates and the vertical black dotted lines represent the corresponding 90% confidence interval for the WGD age estimate, 64–72 Mya. The original raw distribution of dated paralogues is indicated by the circles. The y axis represents the percentage of gene pairs. c, Pruned phylogenetic tree with indication of WGD events (boxes)29. The Cretaceous–Palaeogene (K–Pg) boundary is indicated by an arrow.
We mapped signatures of loss and gain of gene families (Supplementary Note 4.2) onto a phylogenetic tree (Fig. 1a). We also mapped losses and gains of Pfam domains (Supplementary Fig. 4.4, Supplementary Data 6). While many genes are shared between Zostera and Spriodela, clearly some losses and gains are unique to Zostera in relation to its marine environment, the alismatid lineage having set the stage for the subsequent freshwater–marine transition. Those unique to Z. marina include the absence of all the genes involved in stomatal differentiation (Fig. 3a, Extended Data Table 1 and Supplementary Note 5.1) and the disappearance of genes comprising entire pathways encoding volatiles synthesis and sensing (Supplementary Note 6.1), such as those for ethylene15 (Fig. 3b, Extended Data Table 2). Terpenoid genes are also drastically reduced to two (Fig. 3c), as compared with four in Spirodela, 50 in Oryza and > 100 in Eucalyptus, thus precluding synthesis of secondary volatile terpenes (Supplementary Fig. 6.2). Only aromatic acid decarboxylases (AAAD) genes were expanded (Supplementary Fig. 6.3) and these form a clade distinct from Spirodela. The loss of volatiles is also consistent with the loss of stomata, through which they are emitted for airborne communication and plant defence. The repertoire of defence-related genes such as the six groups of NBS_LRR resistance genes (Supplementary Note 6.2) is also reduced to 44 (89 in Spirodela and 100–300 in other plants), which may be linked to a lower probability of infection of Z. marina due to the absence of stomata, which are a main entry point for pests and pathogens in terrestrial plants.
Figure 3: Reconstruction of metabolic (or gene) pathways involved in the production of stomata, ethylene, terpene and pollen in Z. marina.
a, Stomata differentiation from meristemoid mother cells (MMC) to guard mother cell (GMC) to guard cells. b, Ethylene synthesis and signalling up to EIN2 have disappeared; EIN3 and its downstream targets remain. c, Terpenoid biosynthesis in which the pathways producing volatiles are absent but those essential for primary metabolism remain. MVA, mevalonate; MEP, plastidic methylerythritol phosphate; IPP, isopentenyl pyrophosphate; DMAPP, dimethylallyl pyrophosphate; FPP, farnesyl pyrophosphate; GPP, geranyl diphosphate; GGPP, geranylgeranyl pyrophosphate; CPP, copalyl pyrophosphate; GA, gibberellic acid; PP, diphosphate; ABA, abscisic acid. d, Sporopollenin biosynthesis genes; regulatory genes in the nucleus control downstream processes (arrows) in response to signalling coming from external stimuli through receptors on the plasma membrane. All panels: genes in red are absent; blue are present; the grey line represents the plasma membrane. See Extended Data Tables 1, 2, 3.
Land and aquatic floating plants (Embryophyta) are often exposed to intense ultraviolet (UV) radiation and have developed light sensing protein receptors with protective and signalling functions. In contrast, Z. marina inhabits a light-attenuated, submarine environment where it must cope with shifted spectral composition, characterized by low penetration of UV-B, red and far-red wavelengths16. Accordingly, Z. marina has lost ultraviolet-resistance (UVR8) genes associated with sensing and responding to UV damage (Spirodela has not), as well as phytochromes associated with red/far-red receptors (Supplementary Note 7). Whereas photosystems (PSI and PSII) are similar to those of other plants including Spirodela, members of the light-harvesting complex B (LHCB) family are expanded in number, possibly in combination with non-photochemical quenching (NPQ), thereby enhancing performance at low light (Extended Data Fig. 4).
Seagrasses typically experience full marine seawater (35 g kg−1)17, whereas land plants obtain water with low osmolality (0–2 g kg−1) via the rhizosphere and aquatic plants experience fresh (0–5 g kg−1) to brackish (0.5–20 g kg−1) conditions. Although Z. marina displays a typical repertoire of Na+ and K+ antiporters (Supplementary Note 8, Supplementary Table 8.1), one of six H+-ATPase (AHA) genes (Supplementary Table 8.2, Supplementary Data 7) is strongly expressed in vegetative tissue and encodes a salt-tolerant H+-ATPase. Furthermore, Z. marina possesses three AHA genes (along with Spirodela) in a cluster unique to alismatids (Supplementary Fig. 8.1).
Uniquely, Z. marina has re-evolved new combinations of structural traits related to the cell wall. Synthesis of cutin-cuticular waxes to the outside of the leaf epidermis and suberin–lignin near the plasma membrane (Supplementary Note 9, Supplementary Table 9.1) surround a cell wall matrix of (hemi)celluloses, low-methylated pectin (zosterin) and macroalgal-like sulfated polysaccharides18 (Supplementary Note 10). The reduction in carbohydrate-related genes that modify the fine structure of cell wall hemicelluloses and pectins in Z. marina is not due to loss of pathways, but rather to the large variation within these CAZyme gene families in plants. Available genomes (including Spirodela) lack carbohydrate sulfotransferases and sulfatases, suggesting that land plants have lost these genes as a key adaptation to terrestrial as well as freshwater conditions19,20. In contrast, Z. marina has regained the ability to produce sulfated polysaccharides with an expansion of aryl sulfotransferases (12 genes) homologous to aryl sulfotransferases from land plants (Supplementary Note 10). Sulfation facilitates water and ion retention in the cell wall to cope with desiccation and osmotic stress at low tide and, likewise, low methylation of zosterin correlates with the expanded pectin carbohydrate esterase 8 (CE8) family, increasing the polyanionic character of the cell wall matrix. We speculate that several aryl sulfotransferases have evolved because carbohydrate sulfatases have been shown to be active on artificial aryl compounds such as methylumbelliferyl-sulfate21. Osmotic equilibrium is further achieved in Z. marina by organic osmolytes (mainly sucrose, trehalose and proline) in combination with a small cytoplasm:vacuole volume ratio (10%)22. Given that up to 90% of fixed carbon is stored as sucrose in the rhizomes, sucrose synthase (SuSy) and transport (SUT) genes are expanded while those for starch metabolism are greatly reduced, as expected in ‘marine sugarcane’ (Supplementary Note 7.2, Supplementary Data 8).
The repertoire of redox and other stress-resistance genes (Supplementary Note 8) is typical for angiosperms with the exception of catalase (CAT), which is reduced to a single copy in Z. marina (two in Spirodela). Late embryogenesis abundant (LEA) and dehydrins are clearly under-represented in both Zostera and Spriodela relative to other genomes. In contrast, Zostera possesses an unusual complement of metallothioneins. Aside from their role as chelators, metallothioneins may be involved in stress resistance; one of these, MT2L, is among the most highly constitutively expressed genes in Z. marina (Extended Data Fig. 5, Supplementary Note 8.2).
Sexual reproduction of Z. marina takes place underwater, involving completely submerged male and female flowers, and a unique exine-less, filiform pollen that winds around the bifurcate stigmas in a purely abiotic pollination process23. Note that freshwater alismatids (and also Spirodela)24 possess pollen with an exine layer. Exine-less pollen25 is characteristic of all seagrasses except Enhalus acoroides (which is surface pollinated). Ten genes specifically involved in biosynthesis and modification of the pollen exine coat are missing; all other genes involved in the development of viable pollen remain intact (Fig. 3d, Extended Data Table 3, Supplementary Note 11.1). Finally, MADS-box gene transcription factors are also highly reduced to 50 in Z. marina, which is most likely related to its highly reduced flowers (also a feature of Spirodela) that lack the first two whorls of specialized floral leaves, calyx and corolla (Supplementary Note 11.2, Supplementary Table 11.2).
An increasing proportion of the world population inhabits the coastal zone. This impinges multiple pressures on ecosystems including seagrass beds26,27, which in turn compromises the ecosystem services they may provide, including provisioning of harvestable fish and invertebrates, nutrient retention, carbon sequestration and erosion control. In the context of seagrass conservation, elucidating the genomic basis of Z. marina’s complex adaptations to ocean waters (Extended Data Fig. 6) will also inform the development of molecular indicators of their physiological status28, as these unique ecosystems rank, unfortunately, among the most threatened on Earth26,27.
No statistical methods were used to predetermine sample size. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.
A single genotype/clone of Zostera marina (referred to as the ‘Finnish clone’) was harvested on 26 August 2010 at 2 m depth at Fårö Island (latitude 59° 55.234′ N longitude 21° 47.766′ E) located in the northern Baltic Sea, Finland. Plant material was transported to the lab in seawater, cleaned and further processed. Care was taken to use leaf-meristem tissue harvested from the inner layer of basal shoots to minimize bacterial/diatom contamination. Tissues were immediately frozen in LN2 and stored at −80 °C for later DNA and RNA extraction. Monoclonality was verified by genotyping 40 ramets of the mega-clone with six highly polymorphic, microsatellite loci30. There was no evidence for polyploidy25,31,32 (Z. marina is 2n = 12) or somatic mutations33 as assessed by multiple peaks in the microsatellite chromatograms. Tissue was subsequently sent on dry ice to Amplicon Express for HMW DNA extraction using a CTAB isolation method modified by R. Meilan (unpublished) but available from him (rmeilan@purdue.edu), based on the original method34. Following QC according to JGI guidelines, the DNA was shipped to JGI for library and sequencing preparation.
One 35-Kb, fosmid library was generated for end sequencing. The fosmid ends were sequenced with standard Sanger sequencing protocols at the HudsonAlpha Institute for a total of 194,303 Sanger reads (0.29× coverage). Illumina libraries (two fragment libraries (6.62 Gb), one 2-Kb JGI mate-pair library (3.57 Gb), one 4-Kb JGI mate-pair library (3.41 Gb) and two 8-Kb JGI mate-pair libraries (11.94 Gb)) were sequenced with Illumina MiSeq/HiSeq genetic analysers at the Department of Energy’s Joint Genome Institute (JGI), using standard protocols. A total of 25.55 Gb of Illumina and 0.14 Gb of Sanger sequence was obtained representing 47.7× genomic coverage. Prior to assembly, all reads were screened against mitochondria, chloroplast, and Illumina controls. Reads composed of > 95% simple sequence repeats were removed. For the Illumina, paired-end libraries (2 × 250), reads <75 bp were discarded, for the 2 × 150 libraries, reads <50 bp were discarded after trimming for adaptor and quality (q < 20). An additional deduplication step was performed on the mate pairs that identified and retained only one copy of each PCR duplicate. A total of 212,101,273 reads (Supplementary Table 2.1) was assembled using our modified version of Arachne v. 20071016 (ref. 35). Subsequent directed Arachne modules were applied to collapse adjacent heterozygous contigs. The entire assembly was then run through another Arachne process starting at Stage 6 Rebuilder. This produced 15,747 scaffold sequences (30,723 contigs), with a scaffold L50 of 409.5 Kb, 613 scaffolds larger than 100 Kb, and a total genome size of 237.5 Mb (Supplementary Table 2.2).
Scaffolds were screened against bacterial proteins, organelle sequences, GenBank NR (nr_prot) and RefSeq protein databases, and removed if found to be a contaminant. Scaffolds consisting of prokaryotes, chloroplast, mitochondria and unanchored rDNA were removed. We also assembled the chloroplast and partial mitochondrial genomes (Supplementary Notes 2.2 and 2.3, Supplementary Fig. 2.1). Additionally, short (<1 Kb) scaffolds or scaffolds containing highly repetitive sequence ( > 95% 24-mers found more than four times in large scaffolds) or alternative haplotypes were also removed. Following repeat analysis and gene prediction, all scaffolds were subjected to a filtering process (based on NCBI nr_prot + NCBI taxonomy database) to eliminate remaining bacterial (and other) contaminants (Supplementary Table 2.3).
Assembly validation was performed using a set of 12 fully sequenced fosmid clones. In 4 of the 12 fosmid clones, full-length alignments were not found due to fragmentation in the region of the fosmid clone. In five of the remaining eight fosmid clones, the alignments were of high quality (<0.05% bp error). The overall base pair error rate (including marked gap bases) in the fosmid clones that aligned to full length was 0.28% (714 discrepant base pairs out of 253,332 bp). Supplementary Table 2.4 shows the individual fosmid clones and their contribution to the overall error rate. Note that two fosmid clones (16248, 16249) contributed nearly 81% of the discrepant bases. This probably occurred in polymorphic regions of the genome where the haplotype in the fosmid did not match the haplotype in the reference. There are several indels of various sizes in the clone and assembly, typical of a region of degraded transposons. Further quality analysis indicated that 90% of the set of eukaryotic core genes (CEGMA) were present and 98% were partially represented, suggesting near completeness of the euchromatin component.
Two complementary approaches were used to identify repetitive DNA sequences in the Z. marina genome. With respect to masking repeats before gene prediction analysis, a de novo repeat identification was carried out with RepeatModeller (v. open-1.0.7; http://www.RepeatMasker.org)36 to identify repeat boundaries and build consensus models from which potential over represented, non-transposable element, protein-coding genes were removed. RepeatMasker (v. open-4.0.0, WUBlast) was used in combination with this custom repeat library to mask the assembly and prepare it for gene prediction with EuGene.
Furthermore, in order to perform a qualitative and quantitative analysis of repeats with greater resolution37 the genome assembly was processed for de novo repeat detection using the TEdenovo pipeline from the REPET package v. 2.2 (ref. 38); parameters were set to consider repeats with at least five copies. The consensus sequences generated by TEdenovo were then used as probes for whole genome annotation by the TEannot39 pipeline from the REPET package v. 2.2. The consensus repeat sequences were classified using Pastec40. Comparing the genomic positions of transposable elements (TE) to those of exons from the set of predicted genes enabled us to identify that 909 gene predictions most likely represent TEs and these were filtered from the gene set. The REPET package v. 2.2 was also used to annotate repetitive elements in the Spirodela polyrhiza genome assembly with the same parameters as for Z. marina. See Supplementary Fig. 3.1.
Leaf, root and flower tissues were separately frozen in liquid nitrogen immediately following harvest from either ambient (field collected) or experimental (mesocosm) conditions (Supplementary Note 3.2). Overall, we obtained between nine and 20 million high-quality reads from each of the flower-leaf-root replicate libraries; and for the Finnish clone library, 148.5 million high quality reads were retrieved (Supplementary Table 3.3).
The de novo assembly protocol was adapted from ref. 41. We pooled replicates of each tissue together except for the two leaf tissue libraries, which were kept separate (Supplementary Table 3.4) and performed de novo transcriptome assembly for each tissue using Trinity41(v. 2014-07-17) with digital normalization option ON to normalize input read coverage. Frame shift errors and insertion/deletion errors in the assembled transcripts were corrected by FrameDP42. Because a de novo assembly still generates many spurious transcripts, we used the transcript expression value to remove low quality contigs. We used the RSEM pipeline43 to obtain the contig expression values and removed contigs with FPKM (fragments per kilobase of transcript per million fragments mapped) value <1 and IsoPct (percentage of expression for a given transcript compared with all expression from that Trinity component) < 1. In total, we obtained between 39,000 and 53,000 assembled contigs from each library, and 52,000 contigs from the Finnish clone library (Supplementary Table 3.4). Prior to mapping the genome sequence and the predicted genes, we used the CD-HIT44 program (v. 4.6.1) to collapse redundant contigs, which resulted in 79,134 low redundant transcript contigs.
Genomic precursors of known miRNAs were mapped on the Z. marina genome following the procedure described in ref. 47 for the maize genome. miRNA entries from the miRBase database (release 21, 2014) were aligned to the chromosomes of the Z. marina genome. Up to three mismatches were allowed in the alignment, using SeqMap48. In parallel, novel potential DCL1/AGO1-dependent miRNAs were enriched by selecting 5′-U 20–22 nt small RNAs from three different sequenced libraries from Z. marina described in ref. 12. A subset of these small RNAs with abundance ≥10 TPM (transcripts per million) was retained and aligned to the genome with no mismatches. From every locus, we extracted two ~200-nt regions surrounding each aligned miRNA or candidate (from −30 to +160 and from −160 to +30 nucleotides relative to the putative miRNA start or end coordinate, respectively). Minimum energy RNA secondary structures were predicted for each region using the RNAfold program of the Vienna RNA 1.8.5 package (http://www.tbi.univie.ac.at/~ivo/RNA/) using default settings.
In addition, small RNAs from the three sequenced libraries were mapped on these regions, allowing no mismatches, in order to pre-select putative miRNA loci that showed evidence of expression in the three plant tissues analysed. We evaluated RNA structure and small RNA alignment in all the regions based on: (1) dominance of plus-stranded small RNAs; (2) position of the most abundant small RNAs relative to the predicted miRNA coordinates; (3) prevalence of 20–22 nt small RNAs in the predicted miRNA locus; (4) position of the putative miRNA with the stem-loop structure; and (5) absence of oversize (≥3 nt) bulges in the miRNA/miRNA* alignment. After reduction of overlapping loci to a non-redundant set and removal of stem-loop structures with the wrong orientation compared to miRNAs registered in miRBase, we manually inspected the remaining loci to further evaluate them according to the miRNA annotation criteria proposed by ref. 49. Stringency was relaxed when small RNA expression data strongly indicated the presence of miRNA loci that did not meet the whole set of criteria. Novel miRNA precursors overlapping with TEs or other repetitive elements were filtered out.
Potential miRNA targets were identified in silico using the generic small RNA-transcriptome aligner GSTAr from the CleaveLand package (v. 4)50. Predicted targets were accepted with an Allen score <4 or a MFE (minimum free energy) ratio ≥ 7.5. (Supplementary Note 3.4).
Training of the gene prediction programs started with the collection of high quality ESTs. EST information was used, for example, to train the splice predictor SpliceMachine51. Detection of conserved splice sites was further investigated by RNA-seq splice junctions (count > 10) to construct a WAM model in EuGene (v. 4.1)52. Coding-potential was modelled with an interpolated Markov Model (IMM) constructed from the BLASTX alignments of proteins from the PLAZA v. 2.5 database53. An additional protein ‘monocot’ Markov Model was built based on the protein sequences from Brachypodium, maize and sorghum. Starting from EST and protein alignments, a set of 215 gene models was manually constructed and curated using the genome browser GenomeView54. The 215 models were then used as a training set for EuGene in order to optimize the different splice site and coding-potential models, as well as the weights for the extrinsic EST and homology evidence. An overall fitness score of 80.1% was achieved, which is high enough to obtain reliable results without overfitting. GeneMark55 and Augustus56 were separately trained (using the same input data as EuGene) and their predictions were integrated with EuGene using a custom script to evaluate the best gene structure at each locus. All gene models were automatically screened to highlight possible erroneous structures (for example, in-frame stop codons, deviating splice junctions) and manually curated. Transfer-RNA gene models were predicted by tRNAscan-SE (v. 1.31)57 and their structures were verified with Infernal (v. 1.1rc1, rfam11 covariant model database)58. For each gene, UTRs were assigned by identifying a set of ESTs and RNA-seq assemblies that uniquely overlapped with it. We subsequently selected the longest mapped transcript on either end of the predicted coding sequence and designated the section outside the coding sequence as the UTR. Finally, all genes were uploaded to the ORCAE platform (http://bioinformatics.psb.ugent.be/orcae)59, enabling all members of the consortium to refine and curate the gene model and assign gene function. A list of protein domains, as well as the derived Gene Ontology (GO) terms and KEGG pathway identifiers were generated using an InterProScan (v. 5.2.45)60 analysis and available in ORCAE. More specifically, gene functional descriptions were added either manually by consortium expert scientists or automatically through sequence homology searches. The automated method relies on the EC (Enzyme Commission) number reported by InterProScan to retrieve the enzyme name with BLASTP search against UniProtKB/Swiss-Prot61 to filter out hits that are below 60% identity and 70% query/hit coverage. Although such high stringency on per cent identity and sequence coverage reduced the available number of functional descriptions, it reduced the false-positive prediction rate, as desired here.
KS-based age distributions were constructed as previously described62. In brief, the KS values between genes were obtained through maximum likelihood estimation using the CODEML program63 of the PAML package (v. 4.4c)64. Gene families for which KS estimates between members did not exceed a value of 5 were subdivided into subfamilies. For each duplicated gene in the resulting phylogenetic gene tree, obtained by PhyML65, all m KS estimates between the two child clades were added to the KS distribution with a weight 1/m (where m is the number of KS estimates for a duplication event), so that the weights of all KS estimates for a single duplication event summed to one. Mixture modelling was used to confirm a WGD signature in the KS distribution (Fig. 2 and Supplementary Fig. 4.1), for which all duplicates with KS values ≤0.1 were excluded to avoid the incorporation of allelic and/or splice variants, while all duplicates with KS values > 2.0 were removed because KS saturation and stochasticity can mislead mixture modelling above this range62. For further details see Supplementary Note 4.1.
Absolute dating of the identified WGD event was performed as described previously13,29. In brief, paralogueous gene pairs located in duplicated segments (anchors) and duplicated pairs lying under the WGD peak (peak-based duplicates) were collected for phylogenetic dating. Anchors, assumed to be corresponding to the most recent WGD, were detected using i-ADHoRe 3.0 (refs 66,67). Only a low number of duplicated segments and hence anchors could be identified, most likely because of the fragmented assembly of Z. marina. However, the identified anchors did confirm the presence of a broad WGD peak between a KS of 0.8 and 1.6 (data not shown). For each WGD paralogueous pair, an orthogroup was created that included the two paralogues plus several orthologues from other plant species as identified by InParanoid (v. 4.1)68 using a broad taxonomic sampling: one representative orthologue from the order Cucurbitales, two from the Rosales, two from the Fabales, two from the Malpighiales, two from the Brassicales, one from the Malvales, one from the Solanales, two from the Poales, one orthologue from Musa acuminata69 (Zingiberales), and one orthologue from Spirodela polyrhiza11 (Alismatales). In total, about 180 orthogroups from anchor pair duplicates and peak-based duplicates were collected. The node joining the two Z. marina WGD paralogues was then dated using the BEAST v. 1.7 package70 under an uncorrelated relaxed clock model and a LG+G (four rate categories) evolutionary model. A starting tree with branch lengths satisfying all fossil prior constraints was created according to the consensus APGIII phylogeny71. Fossil calibrations were implemented using log-normal calibration priors on the following nodes: the node uniting the Malvidae based on the fossil Dressiantha bicarpellata72 with prior offset = 82.8, mean = 3.8528, and s.d. = 0.5)(ref. 73), the node uniting the Fabidae based on the fossil Paleoclusia chevalieri74 with prior offset = 82.8, mean = 3.9314, and s.d. = 0.5(ref. 75), the node uniting the Alismatales (including Z. marina and Spirodela polyrhiza) with the other monocots based on the oldest fossil monocot pollen, Liliacidites76,77 from the Trent’s Reach locality, with prior offset = 125, mean = 2.0418, and s.d. = 0.5 (refs 14,78) and the root with prior offset = 124, mean = 4.0786, and s.d. = 0.5 (ref. 79). The offsets of these calibrations represent hard minimum boundaries, while their means represent locations for their respective peak mass probabilities in accordance with some of the most recent and taxonomically complete dating studies available for these specific clades14,80. A run without data was performed to ensure proper placement of the marginal calibration prior distributions81. The Markov chain Monte Carlo (MCMC) for each orthogroup was run for 106 generations, sampling every 1,000 generations resulting in a sample size of 104. The resulting trace files of all orthogroups were evaluated manually using Tracer v. 1.570 with a burn-in of 1,000 samples to ensure proper convergence (minimum ESS for all statistics at least 200). In total, 169 orthogroups were accepted and all age estimates for the node uniting the WGD paralogous pairs were then grouped into one absolute age distribution (Fig. 2, too few anchors were available to evaluate them separately from the peak-based duplicates), for which kernel density estimation (KDE) and a bootstrapping procedure were used to find the peak consensus WGD age estimate and its 90% confidence interval boundaries, respectively.
Intra- and inter-genomic co-linearity was investigated (Supplementary Tables 4.2 and 4.3) using MCScanX82 based on a BLASTP search of all genomic protein coding genes with an E-value cut-off of e−10. Only one large duplicated segment was detected, which was most likely due to the fragmented assembly of Z. marina; only 27 scaffolds had a size larger than 1 Mb, accounting for only 23.4% of all protein-coding genes. We therefore additionally used i-ADHoRe (v. 3.0)66 to investigate genomic co-linearity by including all possible scaffolds.
Protein sets were collected for 14 species: Z. marina (ORCAE v. 2.1), Arabidopsis thaliana (TAIR10), Thellungiella parvula (http://thellungiella.org) Populus trichocarpa (Phytozome v. 9.0), Vitis vinifera (Phytozome v. 9.0), Amborella trichopoda (http://amborella.huck.psu.edu), Oryza sativa japonica (Phytozome v. 9.0), Zea mays (Phytozome v. 9.0), Brachypodium distachyon (Phytozome v. 9.0), Spirodela polyrhiza (http://mocklerlab.org), Selaginella moellendorffii (Phytozome v. 9.0), Physcomitrella patens (Phytozome v. 9.0), Chlamydomonas reinhardtii (Phytozome v. 9.0), and Ostreococcus lucimarinus (ORCAE v. 6/3/2013). These species were selected in order to provide a phylogenetic representation traversing green algae, basal plants, monocots, and dicots. Following an ‘all-vresus-all’ TimeLogic Decypher Tera-BLASTP (Active Motif Inc.; e-value threshold 1 × e−3, max hits 500) comparison, OrthoMCL (v. 2.0; mcl inflation factor 3.0)83 was used to delineate gene families. Confidence in establishing gene losses in Zostera was enhanced by using a combination of reciprocal blast, TblastN, re-annotation of Spirodela (and other monocot genes), and careful phylogenetic analysis. OrthoMCL results and related protein resources are available in the ORCAE download section.
To further understand gene family expansion or contraction in Z. marina in comparison with other sequenced genomes, gene family sizes were calculated for all gene families (excluding orphans and species-specific families) (Supplementary Note 4.2). The number of genes per species for each family was transformed into a matrix of z-scores in order to centre and normalize the data. The first 100 families with the largest gene family size in Z. marina were selected. The z-score profile was hierarchically clustered (complete linkage clustering) using Pearson correlation as a distance measure. The functional annotation of each family was predicted based on sequence similarity to entries in the InterProScan and Pfam protein domain database where more than 30% of proteins in the family share the same protein domain. The phylogenetic profile and phylogenetic tree topology provided at PLAZA84 were used to reconstruct the most parsimonious series of gene gain and loss events. The Dollop program from the PHYLIP package85 was used to determine the minimum gene set at ancestral nodes of the phylogenetic tree. The Dollop program is based on the Dollo parsimony principle, which assumes that novel gene families arise exactly once during evolution but can be lost independently in different phylogenetic lineages.
A dedicated search for orthologues/homologues was performed for genes and proteins involved in stomata differentiation (Supplementary Note 5.1), volatile biosynthesis and sensing with focus on ethylene and terpenes (Supplementary Note 6.1), as well as genes involved in male flower specification and pollen differentiation (Supplementary Note 11.1). To this end, queries were chosen from documented genes involved in these pathways (usually from Arabidopsis but occasionally from Oryza, Zea and tomato). Next, the search for homologues in Zostera marina, Spirodela polyrhiza, Oryza sativa japonica and Arabidopsis thaliana (when not used as a query) was performed using BLASTP. To avoid missing or poorly annotated genes a TBLASTN search was conducted using the above queries against the Zostera marina and Spirodela polyrhiza genomes. Putative orthologues were identified based on reciprocal BLASTP searches with Arabidopsis (or the other queries). Owing to species-specific duplications, this may produce some paralogous genes to appear orthologous to the query, or vice versa (see Extended Data Tables 1,2,3). To further confirm correct orthology assignments, phylogenetic trees were built using a broader sampling of protein sequences from both the query species and the three target species. Ambiguously aligned sequences (especially due to indels) were checked manually and corrected or removed.
Raw reads, the assembled genome sequence and annotation are accessible from NCBI under BioProject number PRJNA41721 with GenBank accession number LFYR00000000. The accession number for the Zostera marina Finnish Clone is BioSample SAMN00991190. Fosmid end sequence: GSS KG963492– KG999999; KO000001– KO144970, whole-genome shotgun data: SRA020075 and RNA-seq: GEO GSE67579. Further information on the Zostera marina project is available via the Online Resource for Community Annotation Eukaryotes (ORCA) at http://bioinformatics.psb.ugent.be/orcae/.
Genome sequencing, assembly and automated annotation were conducted by the US Department of Energy (DOE) Joint Genome Institute, Walnut Creek California, USA and supported by the Office of Science of the US DOE, Community Sequencing Program award (2009) contract No. DE-AC02-05CH11231 to J.L.O. Further bioinformatics and annotation was supported in part by the Ghent University Multidisciplinary Research Partnership ‘Bioinformatics: From Nucleotides to Networks’ to Y.V.d.P. Y.V.d.P also acknowledges support from the European Union Seventh Framework Programme (FP7/2007-2013) under European Research Council Advanced Grant Agreement 322739–DOUBLE-UP. RNA-seq (Finnish clone genotype) was funded by the Marine Benthic Ecology and Evolution (MarBEE) group, within the former Centre for Ecological and Evolutionary Studies (now Groningen Institute for Evolutionary Life Sciences), University of Groningen to J.L.O. RNA-seq (flower tissues) was funded by the Excellence Cluster, Future Ocean, Kiel to T.B.H.R. Participation of G.P. and E.D. was supported by the MIUR Italian Flagship project RITMARE (NRP 2011-2013). G.A.P. was supported by FCT-EXCL/AAG-GLO/0661/2012. We thank I. D. Gromicho, KAUST, for his artistry in the production of Extended Data Fig. 6. This work also benefited from discussions within the ESSEM COST action ES0906, “Seagrass productivity from genes to ecosystem management” (2009-2014), J.L.O., G.P. and T.B.H.R; and the Linnaeus Centre for Marine Evolutionary Biology (CEMEB)-Tjärnö, Gothenburg University, J.L.O. and M.T. J.L.O. especially thanks K. Johannesson (CeMEB-Tjärnö), C. Boyen (SBR-Roscoff), R. Reinhardt (MPI-Cologne) and E. Serrão (CCMAR-Faro) for their ongoing encouragement, and the more than 70 colleagues who submitted letters of support for the original proposal to the JGI-Community Sequencing Program.
HudsonAlpha Institute for Biotechnology, 601 Genome Way NW, Huntsville, Alabama 35806, USA.
Marine Ecology Group, Nord University, Postbox 1490, 8049 Bodø, Norway.
Amplicon Express, 2345 NE Hopkins Ct., Pullman, Washington 99163, USA.
J.L.O., T.B.H.R., G.P. and Y.V.d.P. are the lead investigators and contributed equally to the work. J.S., J.W.J., J.G., Y.V.d.P., B.V. and Y.-C.L. coordinated the bioinformatics activities surrounding assembly, quality control, set-up and maintenance of Z. marina on the ORCAE site and deposition of the Z. marina genome resource. T.B.H.R and T.B. generated and analysed RNA-seq libraries from flowers, rhizome, roots. J.L.O., Y.-C.L. and A.J. generated and analysed RNA-seq libraries from the genome genotype and temperature stress experiments. C.B., W.T.S. and J.L.O. contributed to biological sample collection, preparation and quality control prior to DNA extraction. A.M. performed the HMW DNA extraction and quality control from the genome genotype/clone. M.A., J.G., H.T. and M.C. contributed to WGS libraries and sequencing, (fosmid)-cloning and quality control. J.G. coordinated the sequencing of FES, quality control projects. Analysis of architectural features of the genome and annotation of specific gene families, including the written contributions to the main paper and Supplementary Information sections, were performed by the following co-authors: J.W.J., the chromosome assembly analysis; B.V. and Y.-C.L., gene family clustering and comparative phylogenomics; A.R.K. and E.B.B., Pfam domains; E.D.P. and P.J.G., miRNA; R.L., K.V. and Y.V.d.P., whole-genome duplication; F.M., Y.-C.L. and Y.V.d.P., transposable elements; B.V., co-linearity and synteny comparisons; M.T., organellar genomes; P.R., stomata gene family; G.M., cell wall polysaccharides and sulfotransferases; T.T., fatty acid metabolism and its relationship to cell walls and ion homeostasis; P.R., volatiles (ethylene, terpenes); P.R., J.B. and T.B.H.R., metallothioniens; P.R., G.A.P. and C.L., osmoregulation/ion homeostasis/stress-related genes; S.D. and E.D., photosynthetic/ light-sensing genes; G.M., CAZymes; T.B., T.B.H.R. and P.R., plant defence-related; T.B. assembly and analysis of MADS box genes (flowering); P.R.; Y.V.d.P. and Y.-C.L., pollen-related and self-incompatibility genes; F.M., SLR-1gene and core eukaryotic genes analysis (CEGMA). J.L.O., Y.V.d.P., T.B.H.R., C.M.D., Y.-C.L. and P.R. wrote and edited the main manuscript (including the Methods and Extended Data), and organized and further edited the individual contributions (as listed above) for the Supplementary Information sections. J.L.O. and Y.V.d.P. provided the overall evolutionary context and T.B.H.R., G.P. and C.M.D. provided the ecological and societal context. All authors read and commented on the manuscript.
Correspondence to Jeanine L. Olsen or Yves Van de Peer.
Number of genes expressed in five tissues of Z. marina.
Circos plot of the ten largest scaffolds of Z. marina.
Potential impact of transposable elements (TEs) on Z. marina evolution.
Unrooted maximum likelihood tree of genes encoding light-harvesting complex A (LHCA) and LHCB proteins of Z.marina, Spirodela polyrhiza and Arabidopsis thaliana.
Alignment of metallothionein (MT) and half-metallothionein (HMT) genes in Z. marina as compared with other plants.
Conceptual summary of physiological and structural adaptations made by Z. marina in its return to the sea.
This file contains Supplementary Text, Supplementary Figures and Supplementary Tables – see contents page for details.
This zipped file contains Supplementary Datasets 1-8.
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported licence. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons licence, users will need to obtain permission from the licence holder to reproduce the material. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-sa/3.0/.
The seagrass Zostera marina, or eelgrass, is widely distributed throughout the Northern Hemisphere. It is therefore of considerable ecological importance but — as with other seagrasses — its coastal habitats are among the world's most threatened ecosystems. Jeanine Olsen and colleagues report the whole-genome sequence of Zostera. Their analyses provide insights into the evolutionary changes associated with the 'back to the sea' reverse evolutionary trajectory that has occurred in this angiosperm lineage, including the loss of the entire repertoire of stomatal genes, and the presence of sulfated cell-wall polysaccharides that are more macro-algal-like than plant-like.

References: v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v. 
 v.