Transcriptome characterisation and simple sequence repeat marker discovery in the seagrass Posidonia oceanica

Posidonia oceanica is an endemic seagrass in the Mediterranean Sea, where it provides important ecosystem services and sustains a rich and diverse ecosystem. P. oceanica meadows extend from the surface to 40 meters depth. With the aim of boosting research in this iconic species, we generated a comprehensive RNA-Seq data set for P. oceanica by sequencing specimens collected at two depths and two times during the day. With this approach we attempted to capture the transcriptional diversity associated with change in light and other depth-related environmental factors. Using this extensive data set we generated gene predictions and identified an extensive catalogue of potential Simple Sequence Repeats (SSR) markers. The data generated here will open new avenues for the analysis of population genetic features and functional variation in P. oceanica. In total, 79,235 contigs were obtained by the assembly of 70,453,120 paired end reads. 43,711 contigs were successfully annotated. A total of 17,436 SSR were identified within 13,912 contigs.


Background & Summary
Posidonia oceanica is an endemic seagrass species in the Mediterranean Sea where it forms extensive meadows along the coastline, from the surface to 40 m depth 1 . P. oceanica provides important ecosystem services including sediment stabilization, contribution to CO 2 absorption, sheltering and food sources for a diverse associated community, which also includes commercially important fish species 2 . The total mapped extent of P. oceanica is estimated about 1.2E+6 ha 3 with an estimated economic value of 1.7E+6 € ha y − 1 (ref. 2). P. oceanica meadows are declining throughout the Mediterranean basin, due to direct and indirect effect of human activities 3 , with important consequences on the stability of the coastline and on the whole set of ecosystem services provided by meadows. Plants can respond to changes adopting meadow-scale adaptive strategies or can exhibit organismal-scale responses, adjusting fitness and performance. Understanding how human-driven environmental changes impact marine plants and how plants respond to changes is critical to design appropriate conservation strategies. However, our understanding of the consequences of the human impact on P. oceanica meadows is hampered by the absence of specific molecular tools. To date, a quite extensive literature exists on the genetic diversity and structure of the species, using a limited set of SSR (Simple Sequence Repeats) markers e.g., refs 4-6. Nevertheless, the limited number of markers and the neutral behavior of most of them, did not allow to detect clear signature of selection along environmental and geographic gradients. Only a limited number of functional studies has been conducted on P. oceanica, mostly looking at the expression of selected target genes e.g., refs 7-10. A first set of annotated genes was developed for P. oceanica, reconstructed from samples collected at different depths, and a SSH (Suppression Subtractive Hybridization) library was built in order to detect differentially expressed genes 7,8,11 . These resources constitute an important first step towards the characterization of the transcriptome of the species, but the limited number of sequences is not at all representative of the whole transcriptome of the species. Here we make a step forward, building a complete transcriptome that can be interrogated to investigate the effect of human-driven environmental change on the endemic seagrass species P. oceanica, also taking advantage of the recent release of the first seagrass genome from the temperate species Zostera marina 12 .
P. oceanica is a clonal grass species, with lignified rhizomes and a bundle of erect strip-like leaves. Sexual reproduction is not regular and the flower is hermaphroditic. Shallow and deep portions of the same meadow are usually genetically distinct (e.g., ref. 13). We generated a de novo transcriptome assembly by sequencing specimens collected at two different depths and at different times of the day to capture representative transcriptional changes associated with photoperiod and habitat distribution in our species. The assembled contigs were annotated, generating the first gene catalogues for P. oceanica. In addition we interrogate the newly assembled transcriptome to identify Simple Sequence Repeats. Next generation sequencing (NGS) technologies recently become the methods of choice for gene expression analysis due to their unprecedented level of sensitivity and high-throughput nature 14 . RNA-seq has also been shown to be a very efficient, cost-effective approach for transcriptome analysis of non-model species. In fact, also in absence of the genome, RNA-seq allows construction of complete transcriptome of an organism by de novo assembly 15 .
The catalogue of annotated P. oceanica genes together with the SSR provides the seagrass community with advanced tools to link the effect of human-induced environmental change and functional responses in P. oceanica, and molecular markers useful in population genetic studies. These resources will be of critical importance to implement the conservation strategies aimed at preserving this iconic species.

Sampling design and RNA extraction
Posidonia oceanica shoots were collected in 2012 from a meadow located in Stareso (Corse, 8°45′E, 42°35'N, Fig. 1). Sampling was performed at two different depths, −5 m (shallow station) and −25 m (deep station) and at two different times of the day (day and night). Sampling was performed in the framework of a comprehensive study aimed to characterize P. oceanica meadows, from the plant to the whole ecosystem 16 , in collaboration with international partners within the EU COST Action ES-0906 (Seagrass productivity: from genes to ecosystem management). Leaf tissue from an adult shoot from each depth and day-time was cleaned from epiphytes, stored in RNA later solution (Life Technologies) overnight at 4°C and then transferred at −20°C. Total RNA was isolated using the Aurum Total RNA mini kit (Bio-Rad) following the manufacturer's instructions. The quality of total RNA was checked using NanoDrop technologies (ND-1000 Spectrophotometer, Peqlab) and the Agilent 2100 Bioanalyzer (RNA Nano Chip, Agilent).

Transcriptome sequencing, assembly and annotation
Libraries were constructed starting from 3.5 μg of total RNA using the TruSeq RNA Sample Preparation Kit. Paired ends libraries (95 bp × 2) were sequenced on an Illumina Genome Analyzer IIx instrument (Illumina, San Diego, California, USA) following manufacturer instructions. Prior to assembly, raw reads were filtered with Prinseq v0.20.4 (ref. 17), in order to remove duplicated sequences and the reads showing a quality score mean below 30 (Qo30) and with Cutadapt v1.9, to remove adapters (Table 1). High quality reads were pooled and used for the de novo assembly using Trinity v2. 1 20) was used to calculate the coverage depth. Contigs were cleaned from redundant sequences using CD-HIT program v4.6 (ref. 21). To assess the completeness of the transcriptome assembly, we used the 458 Core Eukaryotic genes (KOGs) obtained from CEGMA (http://korflab.ucdavis.edu/datasets/ cegma/) to measure how the (KOGs) set is covered by the P. oceanica transcriptome. Analysis was performed using the TransRate software 22 which makes use of the reciprocal best blast hit methods to compare a transcriptome against a reference. The obtained contigs were annotated by blast search (BLASTx) against the NCBI nr protein database or by conducting a local BLASTx search. The databases queried with a local BLASTx search included the Vitis vinifera proteins, Phoenix dactylifera proteins and Elaeis guineensis proteins, downloaded from NCBI GenBank, Zostera marina proteins downloaded from ORCAE 23 Clusters of Orthologous Group (COG) 24 . Information about domain/motifs content of the sequences were retrieved using the InterProScan functionality in Blast2GO v3. 1.3 (ref. 25). To classify the function of contigs, GO assignment was performed using Blast2GO software.

Detection of SSR containing sequences.
SSR were identified using the MISA MicroSatellite identification tool (http://pgrc.ipk-gatersleben.de/ misa/), defining a minimum repeat length criteria of six repeated units for dinucleotides, five repeated units for tri, tetra, penta and hexanucleotides. We also identified compound SSR as the sequences in which two SSR were separated by no more than 10 bases.

Data Records
The project was submitted to NCBI BioProject with BioProject ID: PRJNA315106 (Data Citation 1). The assembled contigs were deposited at DDBJ/ENA/GenBank under the accession GEMD01000000 (Data Citation 2).

Technical Validation
In order to characterize and comprehensively cover the transcriptome of Posidonia oceanica, RNA was isolated from plants living at different depth (shallow and deep) and in different times of the day (day and night). Four libraries were prepared and sequenced using the Illumina sequencing technology. A total of 86,691,134 reads were generated. After filtering, 70,453,120 high quality reads were obtained and used for de novo assembly (Table 1). Reads were assembled using the Trinity software obtaining a total of 165,235 contigs. The reads used to create the assembly were mapped back to the assembly using Botwie and the coverage depth calculated. Contigs showing a coverage less than a threshold of 5X were removed. CD-HIT clustering program was used to remove sequence redundancy by applying a global identity threshold of 90% and retaining the longest possible contigs. After these filtering steps, 79,235 contigs  remained. As shown in Table 2, the contigs lengths ranged from 201 to 16,705 bp, with an average length of 1,278 bp. The contigs size distribution is shown in Fig. 1 (Table 3). These results suggest that the assembly is a reliable representation of the P. oceanica transcriptome. Functional annotation of the P. oceanica transcriptome was based on two levels of sequence similarity namely sequence-based and domain-based alignments 26 . Sequence based alignments were performed by blast search (BLASTx) against the non-redundant protein database (NR) in NCBI and by local BLASTx against the proteins of the four top hit species identified from the BLASTX against the NR protein database. In Fig. 1 the effect of the sequence length on the percentage of significant matches in BLASTx results is reported. The data shown in Fig. 1 make evident that the longer the length of the assembled sequences the greater the number of matches. The E-value distribution in Fig. 2a revealed that    valueo1.0e-50) while the remaining 40.2% show a match with an E-value above 1.0e-50. In Fig. 2b, the BLASTx top hit species distribution of matches with known sequences is reported and indicate that the majority of P. oceanica contigs show the highest homology with Z. marina sequences (24.39%). The others most represented species included E. guineensis (14.21%), P. dactylifera (11.91%) and V. vinifera (8.79%). All the alignments have been carried out setting the E-value thresholds at the value of ≤1e-5. Annotation was completed aligning the contigs against COGs and using the functionality of InterProScan in Blast2GO software. 22,035 assembled contigs matched entries against the InterProScan queried protein databases (Table 4). A local BLASTx was used to align the assembled contigs against the COGs database 21 to predict and classify possible function. 17,421 COG functional annotations were obtained that were grouped into 24 functional categories (Fig. 3). The annotation results are summarized in Supplementary Table S1.
The Fig. 4 shows Venn diagrams illustrating the distribution of similarity search results. In the Venn diagram in Fig. 4a the BLAST results against NR protein, Z. marina proteins, E. guineensis proteins, P. dactylifera proteins and V. vinifera proteins are reported. The total number of annotations obtained     from these databases was 43,325 (54.6% of total final contigs). In Fig. 4b, the annotation results obtained from COGs and InterProScan are shown. The number of contigs which show a hit against these databases was 27,646 (34.8% of total final contigs). In the Venn diagram in Fig. 4c, all the annotation resulting from the databases in Fig. 4a,b are reported. Integrating the results obtained from all the databases queried, a total of 43,682 sequences (55% of assembled contigs) resulted annotated. Blast2GO PRO was used to assign Gene Ontology (GO) term and functionally categorize the assembled P. oceanica contigs. At least one GO is assigned to 25,495 assembled contigs (32.1% of total contigs; Fig. 5).
The annotated transcript sequences here reported represent a significant improvement of the genomic information available for this species.
The 79,235 contigs obtained from the assembly were mined for the identification of contigs containing SSR. The search criteria adopted allowed the identification of 17,436 SSR contained in 13,912 contigs (Supplementary Table S2 and Table 5). The most abundant repeat type was trinucleotides (58%) and dinucleotides (38%). The most abundant repeat motif was AG among dinucleotides (18% of the dinucleotides repeats), GGA among trinucleotides (6.4% of the trinucleotides repeats), AAAG among tetranucleotides (8.3% of the tetranucleotides repeats). Once the SSR here discovered will be tested for polymorphism, this will increase the SSR markers resources for this species and will be of help in population genetics studies. Furthermore, because these SSR were found in the transcribed part of the genome they will greatly help the studies aimed to identify the genes involved in the adaption of P. oceanica in the deep-sea environment.