Genomic analysis of the native European Solanum species, S. dulcamara

De novotranscriptome assembly

Short reads from seventeen different S. dulcamara cDNA libraries that were sequenced using either Roche GS-FLX or Illumina HiSeq2000 sequencing technologies (Table 1; Additional file 1: Table S1) were combined to build de novo a consensus transcriptome using the Trinity package [14]. This resulted in an assembly of 32,157 contigs of more than 500 nts in size, with an average length of 1,346 nts (Table 2; Additional file 1: Table S2). The dataset encompasses 24,193 unigenes, of which 3,885 are clusters with multiple variants (Table 2; Figure 2). These variants are expected to comprise allelic variants, splice variants, nearly-identical paralogs or mis-assemblies. The sequences of all contigs are available at the Sol Genomics Network website (http://ftp.solgenomics.net/unigene_builds/single_species_assemblies/Solanum_dulcamara/).

Library	Sequencing method	# Reads (raw)	# Reads (high quality)	Avg read length (high quality)
Mixed1	GS-FLX titanium, single end, random primed, normalised	799,937	786,113	297
Mixed1	GS-FLX titanium, single end, 3’ primed, normalised	568,762	558,729	312
Leaves1	HiSeq2000, single end, random primed, 100bp, normalised	127,318,285	107,557,660	96
Stem & primoridia1	HiSeq2000, single end, random primed, 50bp	516,292,741	471,046,908	50

¹ See Additional file 1: Table S1 for details.

Total # of contigs	32,157
Total # of unigenes	24,193
# consisting of single sequence	20,308
# consisting of multiple variants	3,885
Total sequence length (nt)	43,277,997
Average contig length (nt)	1346
Minimum contig length (nt)	501
Maximum contig length (nt)	15047
Median contig length (nt)	1043

Functional annotation

Comparison of protein family structure between S. dulcamaraand other plant species

Multi-species transcriptome comparison may be used in order to identify orthologous gene groups, measure changes in the size of protein-coding gene families, study gene family evolution and detect taxonomically restricted sequences (i.e. species-specific or genus-specific sequences).

ORF/protein prediction

To be able to compare protein family structure between S. dulcamara and other plant species we first predicted the ORFs and protein sequences encoded by the S. dulcamara contigs. ESTScan [20] of the 32,157 contigs indicated that 26,696 contigs (~83%) contain putative coding sequences that could be translated into proteins. This is very similar to the percentage of contigs predicted to be protein coding by BLASTx (see above), with the slightly higher percentage of the latter probably explained by the fact that BLASTx better tolerates sequencing errors that result in frame shifts and premature stop codons than ESTScan. In total, 11,760 full-length proteins (i.e. containing a putative start and stop codon) and 14,936 truncated proteins (from partial coding sequences) were identified (Table 3). To confirm the reliability of the ESTScan prediction we carried out BLASTp searches of the predicted proteins against the tomato, potato and Arabidopsis protein complement. About 95% (25,418) of the S. dulcamara proteins had a significant match in at least one of these protein databases. Comparison of the BLASTp results with the BLASTx results of the same contigs revealed that in 99.9% of the cases, the best hit was identical. As a measure of the quality of our assembly, we also compared the size distribution of the subset of S. dulcamara full-length proteins to the length distribution of the proteins encoded in the genomes of tomato (34,727) and potato (39,031), the two Solanum species for which a full genome sequence was published recently [12, 13]. Although the number of S. dulcamara full-length proteins is three to four times smaller than the number of proteins in the tomato and potato genome, protein size in the three datasets shows a similar log-normal distribution (Figure 3). Together, these results support the reliability of the assembly and the predicted protein data set.

Items	# contigs
complete ORF	11,760
5'truncated ORF1	8,608
3' truncated ORF2	3,405
5'and 3' truncated ORF3	2,923
no good ORF4	5,461
Total	32,157

¹ ORFs lacking the initiating ATG codon, but including a termination triplet.

² ORFs including the ATG start codon, but lacking the stop codon.

³ ORFs having neither start nor stop codon.

⁴ Transcripts showing interspersed stop codons.

OrthoMCL clustering

Orthologous gene groups were identified using orthoMCL [21]. The analysis included protein datasets from S. dulcamara, from the related Solanum species tomato and potato, as well as from the two model plant species Arabidopsis and rice. As the input for S. dulcamara we used the partial and full-length proteins predicted by ESTScan. To ensure that each locus was represented only once in the orthologous gene group analysis, only the longest predicted protein from each variant cluster was used. Similarly for the other species, only the longest protein variant encoded by a locus was used. A total of 164,689 protein sequences from the five species were clustered into 23,370 ortholog groups (see Additional file 1: Table S3). A consensus annotation was automatically assigned to each group based on the frequency (greater or equal to 0.33) of the most prevalent InterPro entry list. In case the threshold criterion was not met, the combination of the two most frequent InterPro entry lists was used. In Figure 4, the number of orthologous and putative species-unique gene groups is shown. Of the 19,713 proteins from S. dulcamara, 15,073 were placed in a total of 13,518 gene groups with multiple members and 4,640 were not grouped and defined as species-specific singletons. As expected, a large part of the S. dulcamara gene groups (7,737) contained orthologs from all other species, thus representing genes that are highly conserved in flowering plants. High sequence conservation and high gene expression have been suggested to correlate [22–24], which may explain why the RNAseq-based S. dulcamara transcriptome has a slight bias towards highly conserved gene groups (43%), compared to the transcriptomes of tomato (31%) and potato (32%), which were derived from whole-genome sequencing. In S. dulcamara, as in the other species, many genes were species specific: 17 gene groups and 4,640 singletons.

SSR analysis

In-silico SSR prediction

As a first effort aimed at the development and characterization of EST-based SSR markers in bittersweet, the S. dulcamara contigs were examined for the presence of SSR motifs using MISA (http://pgrc.ipk-gatersleben.de/misa/; Additional file 1: Table S4). A total of 6,029 SSRs were identified (Table 5). The frequency and the number of repeat units for each of the SSR motifs (except the mono-type repeats) are shown in Additional file 1: Table S5. With the exception of mono-type and compound repeats, the average length of the SSRs was ~16 nucleotides (range from 12 to 96 nucleotides). By exploiting the ORF/protein predictions from ESTScan, we also investigated the position of the SSRs along the transcripts (Figure 5). Tri-nucleotide as well as hexa-nucleotide repeats are preferentially located in the coding regions, while the other types are more frequent in UTRs. In particular, di-, tetra- and penta-nucleotide are preferentially located in 5' UTRs, whereas mono-nucleotide SSRs are equally distributed between 5' and 3' UTRs. These results are in agreement with observations in [25]. Using the Primer3 software [26], primer pairs to amplify each SSR were successfully designed for 4,233 transcripts and failed for 1,537 sequences (see Additional file 1: Table S4).

Items	Counts
Total number of identified SSRs:	6,029
Number of SSR containing sequences:	5,156
Number of sequences containing more than 1 SSR (c)	742
Mono-nucleotide (p1)	3,060
Di-nucleotide (p2)	1,163
Tri-nucleotide (p3)	1,733
Tetra-nucleotide (p4)	42
Penta-nucleotide (p5)	9
Hexa-nucleotide (p6)	22

Genetic map construction

S. dulcamara contains a number of traits, such as several pathogen resistances, which may be valuable for agricultural purposes. To speed up gene mapping efforts in this species we generated a genetic map that is anchored to the high-quality genetic/physical maps of tomato. In addition to providing practically useful information on synteny and co-linearity between these two species, such a map also offers insight into genome evolution in the genus Solanum.

Segregation analysis and map construction

A subset of SNPs was identified that was heterozygous in the female parent and homozygous in the male parent (see methods for criteria) and 96 of these were selected in such a way that the tomato orthologs of the corresponding contigs were distributed evenly over the tomato genome (Additional file 1: Tables S2 and S8). Segregation of 90 out of 96 selected SNP could successfully be determined in 94 F₁ individuals, using KASPar assays (http://www.lgcgenomics.com/). In addition, segregation of 108 AFLP and 27 CAPS markers was analysed in the same individuals. Linkage analysis using JoinMap [31] revealed the existence of 12 linkage groups, in agreement with the haploid chromosome number of the species. The resulting linkage groups ranged from 76 to 121 cM in size and harboured 10 to 30 markers each (Figure 6; Additional file 1: Table S9). Synteny and co-linearity with tomato was studied using marker orthology (Figure 7; Additional file 1: Table S8). To better understand chromosome evolution in the Solanaceae, the genetic map was subsequently compared to the integrative maps of tomato (S. lycopersicum), potato (S. tuberosum), eggplant (S. melongena), pepper (Capsicum annuum), tobacco (Nicotiana) and their deduced common ancestors, as presented by Wu and Tanksley [32] (Additional file 2: Figure S1).

Plant material

S. dulcamara material used to generate mRNA samples for RNAseq is described in Additional file 1: Table S1. Material used to test SSRs was provided by Dr Janny Peters (Radboud University Nijmegen, The Netherlands). The segregating population used for map construction (code 05-150) was derived from a cross between accession A54750069-1 and 944750001-2 [10]. All plants were cultivated in standard greenhouse conditions as described in [10], unless indicated otherwise.

RNA extraction and sequencing

Total RNA was isolated using Trizol (for the “Mixed” and “Leaves” libraries, see Table 1) or the Plant RNeasy kit (Invitrogen; for the “Stem & primordia” library, see Table 1) and treated with DNase. In case of the “Mixed” libraries, mRNA was purified and duplex-specific-nuclease-normalized cDNA samples were prepared and sequenced by Eurofins MWG Operon (Ebersberg, Germany) on the Roche GS-FLX platform. In case of the “Leaves” library, mRNA was purified and duplex-specific-nuclease-normalized cDNA samples were prepared and sequenced by Fasteris SA (Geneva, Switzerland). For the “Stem & primordia” library, mRNA was purified and cDNA samples were prepared and sequenced by Fasteris SA without prior normalisation.

De novotranscriptome assembly

Raw read filtering based on quality values and length was performed with the ‘Trim sequences’ algorithm in CLC Genomics Workbench v4.7.1 (CLCBio, Aarhus, Denmark). Default settings were used and low quality sequences (limit=0.05) and sequences no longer than 50 nts were removed. Although the assembler algorithm discarded low coverage k-mers, the raw reads were error corrected in order to speed up the assembly process. Therefore all sequence data except the Roche GS-FLX data was base-error-corrected with decGPU version 1.06 [35]. DecGPU was run with default settings. The decGPU algorithm output consisted of error free reads, fixed reads and discarded reads. For the assembly both error free and fixed reads were used. The decGPU process discarded 66M sequences (11% of total Illumina input sequences). All samples where pooled, both Roche GS-FLX and Illumina sets, and assembled using the de novo transcriptome assembler Trinity version 2011-10-29 (http://trinityrnaseq.sourceforge.net/) [14]. The Trinity assembly was run with a default fixed k-mer length of 25, minimal contig length of 500 bp, minimal k-mer coverage of 2 and a butterfly heap space size of 50GB.

ORF identification and functional annotation

Automated annotation was performed by BLASTp and BLASTx searches (e-value < e^-10) against the S. lycopersicum (version iTAG 2.3 proteins), S. tuberosum (version 3.4), A. thaliana protein complement (version TAIR 10 pep) and the UniProtKB/Swiss-Prot database (release 2012_02). In addition, BLASTn searches (e-value < e^-10) against the nucleotide non-redundant database (GenBank, release March 2011) were carried out. The Blast2GO suite (version 2.5.0) [18] was used to identify InterPro entries that were mapped to GO terms. KAAS [19] was used to assign KO (KEGG orthologs) terms (representative gene data set for eukaryotes + plants) to S. dulcamara transcripts. The BBH (bi-directional best hit) option was used to map KO terms onto KEGG pathways, using the same program.

Identification and annotation of orthologous gene groups

ESTScan [20] was used to predict ORFs in the S. dulcamara transcriptome using the default Arabidopsis thaliana training matrix for peptide prediction. OrthoMCL (version 2.0.2) [21] was used to identify gene family groups among S. dulcamara (19,713), S. lycopersicum (34,727), S. tuberosum (39,031), A. thaliana (27,416), O. sativa (43,802). Enclosed within brackets, is reported the number of proteins used as input data, after removing all but the longest protein sequence in case of splice variants. All the resulting sequences were merged into a single FASTA file and all-versus-all comparisons were performed using BLASTp (e-value < e^-5). For the MCL clustering algorithm we used an inflation value (-I) of 1.5 (OrthoMCL default). Consensus annotation of each gene group was automatically assigned based on of the most frequent InterPro entry list (frequency ≥ 0.33). In case the threshold criterion was not satisfied, the combination of the two most frequent InterPro entry lists was used. In case of Arabidopsis, rice and tomato we exploited the already available nterPro annotations (http://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_functional_descriptions; http://rapdb.dna.affrc.go.jp/download/archive/irgsp1/IRGSP-1.0_representative_2012-05-28.tar.gz; http://ftp.solgenomics.net/tomato_genome/annotation/ITAG2.3_release/ITAG2.3_desc_and_GO.csv). In contrast, because no InterPro annotation is available at http://potatogenomics.plantbiology.msu.edu/index.html, we identified the InterPro protein domains within the potato sequence collection using the Blast2GO suite. The GO term enrichment analysis was performed using TopGO package (v2.8) from Bioconductor (http://www.bioconductor.org/packages/2.12/bioc/html/topGO.html). Fisher's exact test (p-value < 0.01) was used to identify the over-represented GO terms.

SSR identification and analysis

The SSR search tool MISA (MIcroSAtellite identification tool; http://pgrc.ipk-gatersleben.de/misa/) was used to identify and localize single or multiple stretches of microsatellite motifs. Research criteria include a minimum of 10 in case of mononucleotide and a minimum of 4 repetitive units in case of 2-, 3-, 4-, 5-, 6- unit repeats. Primer pairs flanking the microsatellite loci were automatically designed using Perl scripts provided with MISA (http://pgrc.ipk-gatersleben.de/misa/primer3.html) in combination with the software Primer3 (http://www.broadinstitute.org/genome_software/other/primer3.html) [26].

Genomic DNA was isolated from seven individuals using a standard CTAB method [36, 37]. PCR amplifications were performed in 20 μL volume and the reaction mixture contained 16 ng of genomic DNA, 100 μM of each primer, 400 μM of dNTPs each, 0.5 units of DreamTaq Polymerase (Fermentas), 1× DreamTaq DNA polymerase buffer containing 20 mM MgCl₂. PCR conditions consisted of an initial denaturation step at 95°C for 2 min, followed by 40 cycles of 95° C for 20 sec, 55°C for 30 sec, and 72°C for 20 sec and a final step at 72°C for 30 min. Forward primers were labelled with WellRED D2-PA fluorescent dye (Sigma-Aldrich). 1 μL of PCR product and 0.3 μL of GenomeLab DNA Size Standard 400 (Beckman Coulter) were diluted with 30 μL sample loading solution (Beckman Coulter) and electrophoresis was performed on the CEQ 8800XL Genetic Analysis System (Beckman Coulter). SSR locus allele sizes were determined using the Beckman CEQ fragment analysis software.

SNP discovery and KASPar marker design

For SNP discovery, RNAseq sequence reads were first trimmed and mapped to the S. dulcamara transcriptome assembly using CLC Genomics Workbench v4.7.1 (CLCBio, Aarhus, Denmark). For trimming, low quality sequence (limit 0.05) and ambiguous nucleotides (if longer than 2 nucleotides) were removed and reads shorter than 50 nucleotides were discarded. As a template, we selected contigs with an average coverage of more than 20 and less than 10,000 reads. For read mapping, the minimal read coverage was set at 90% and the minimal alignment identity at 90%. Reads that could be mapped to multiple locations with the same score (repeats) were assigned randomly to one of these locations. SNP calling was done using an upgraded version of QualitySNP [38] (H. Nijveen, Wageningen University, The Netherlands, unpublished). To determine valid SNPs, the minimum similarity score per polymorphic site was set at 0.75 and the minimum similarity score of all polymorphic sites at 0.8, INDEL SNPs were marked as low quality and removal of low quality at sequence ends was disabled. For SNP selection, SNPs were taken as heterozygous in A54750069-1 when coverage was at least 10 reads and frequency of each allele was more than 30%. SNPs were taken as homozygous in 944750001-2 when coverage was at least 20 reads and the number of alternative allele reads was 0. Furthermore, the GC content of the 50bp flanking the SNP on either side was between 30 and 70% to aid in primer design. Primers for KASPar assays (http://ftp.solgenomics.net/tomato_genome/annotation/ITAG2.3_release/ITAG2.3_sgn_data.gff3) were designed by KBioscience (Hoddesdon, UK) and assays were performed according to the manufacturer’s protocol on a Fluidigm EP1 system (Fluidigm Corporation, San Francisco, CA). A BLASTp-based analysis was performed to identify probable orthologous pairs as reciprocal best hits between bittersweet and tomato. BLASTp hits with bit scores lower than 100 were filtered out. In addition, BLASTp searches were performed to identify near identical paralogs within the reference species tomato. In case the bit score of a paralogous pair was higher than the score of the corresponding orthologous pair, the corresponding tomato transcript was regarded as having in-paralogs (i.e. paralogs arising from duplication after speciation) and, unless paralogs were located not more than 10 loci apart, discarded from the analysis. The position of tomato genes on the genetic map was estimated through the identification of the closest 3’ and 5’ marker from the tomato Kazusa F2-2000 linkage map available at the SGN ftp site [39]. Physical positions were retrieved from a gff3 file that contains alignments of marker sequences to the tomato pseudo-molecules (http://solgenomics.net/itag/release/2.3/list_files/ITAG2.3_sgn_data.gff3).

Genetic map construction and comparison

An S. dulcamara F₁ population containing 94 individuals was screened with 328 SNP, CAPS and AFLP markers (Additional file 1: Tables S8 and S9). Markers with >75% missing data points, completely linked markers (except for SD2/SD58, used for anchoring to the tomato map) and markers showing a segregation ratio significantly different from 1:1 (χ² test, p<0.005) (except for SD2/SD58, used for anchoring to the tomato map) were excluded prior to analysis. Regarding the KASPar assays, in one case no SNP was detected and in four cases the assay did not work properly. Matrix data was analysed using JoinMap v4.1 [31] with CP (cross-pollination) population type settings. For linkage analyses the regression and maximum-likelihood mapping algorithms were used, with Haldane’s mapping function. Markers with different relative positions on the two generated maps and markers that could not be assigned to a linkage group were rejected. The positions of 225 markers, as obtained by regression mapping, are presented here. Rooting of groups was done with LOD>4. For comparison of the S. dulcamara genetic map with those of other Solanaceae as produced by Wu and Tanksley [32], all markers were mapped to the tomato Kazusa F2-2000 linkage map. Three markers were excluded from the genetic map because they suggested singleton translocations that were not corroborated by neighbouring markers: Agl (S. dulcamara chr1, 65.0 cM, S. lycopersicum chr4, 36.5 cM), SD27 (S. dulcamara chr11, 57.0 cM, S. lycopersicum chr4, 37.4 cM) and SD122 (S. dulcamara chr6, 63.0 cM, S. lycopersicum chr2, 35.2 cM).

Accession numbers

Raw sequence reads obtained from 454 and Illumina sequencing were submitted to the NCBI Short Read Archive (SRA, http://www.ncbi.nlm.nih.gov/sra) under the accessions SRP020226. The S. dulcamara contigs are available from the Sol Genomics Network website (http://ftp.solgenomics.net/unigene_builds/single_species_assemblies/Solanum_dulcamara/) and are included in the SGN BLAST search tool (http://solgenomics.net/tools/blast/index.pl?db_id=203).