Reference genomes and transcriptomes of Nicotiana sylvestris and Nicotiana tomentosiformis

Genome assembly

The N. sylvestris and N. tomentosiformis genomes were sequenced using a whole-genome shotgun sequencing approach. For N. sylvestris, a 94× coverage of 100 bp Illumina HiSeq-2000 reads was used. In total, six libraries were constructed with different insert sizes ranging from 180 bp to 1 kb for paired-end libraries, and from 3 to 4 kb for mate-pair libraries. The numbers of clean reads in each library are summarized in Additional file 1. Similarly, for N. tomentosiformis a 146× coverage of 100 bp Illumina HiSeq-2000 reads was used. In total, seven libraries were constructed with different insert sizes ranging from 140 bp to 1 kb for paired-end libraries, and from 3 to 5 kb for mate-pair libraries. The numbers of clean reads in each library are summarized in Additional file 2.

With the advent of next-generation sequencing, genome size estimations based on k-mer depth distribution of sequenced reads are becoming possible [36, 50–52]. For instance, the recently published potato genome was estimated to be 844 Mb using a 17-mer distribution [36], in good agreement with its 1C size of 856 Mb [4]. Furthermore, the analysis of repetitive content in the 727 Mb potato genome assembly and in bacterial artificial chromosomes and fosmid end sequences indicated that much of the unassembled genome sequences were composed of repeats [36]. In N. sylvestris and N. tomentosiformis the genome sizes were estimated by this method using a 31-mer to be 2.68 Gb and 2.36 Gb, respectively. While the N. sylvestris estimate is in good agreement with the commonly accepted size of its genome based on 1C DNA values, the N. tomentosiformis estimate is about 15% smaller than its commonly accepted size. Estimates using a 17-mer were smaller; 2.59 Gb and 2.22 Gb for N. sylvestris and N. tomentosiformis, respectively. Using the 31-mer depth distribution, we estimated that our assembly represented 82.9% of the 2.68 Gb N. sylvestris genome and 71.6% of the 2.36 Gb N. tomentosiformis genome.

The proportion of contigs that could not be integrated into scaffolds was low; namely, the N. sylvestris assembly contains 59,563 contigs (18 Mb; N50 length of 273 bp) that were not integrated in scaffolds, and the N. tomentosiformis assembly contains 47,741 contigs (17.3 Mb; N50 length of 346 bp) that were not integrated in scaffolds. Using the regions of the Whole Genome Profiling (WGP) physical map of tobacco [53] that are of N. sylvestris or N. tomentosiformis ancestral origin, the assembly scaffolds were superscaffolded and an N50 of 194 kb for N. sylvestris (10,261 contigs and scaffolds in 2,637 superscaffolds) and of 166 kb for N. tomentosiformis (7,463 contigs and scaffolds in 1,989 superscaffolds) were obtained. Superscaffolding was performed using the WGP physical map contigs as templates and positioning the assembled sequences for which an orientation in the superscaffolds could be determined. This approach discards any anchored sequence of unknown orientation as well as any sequence that spans across several WGP contigs, thereby reducing the number of superscaffolded sequences. Furthermore, the superscaffolding introduced additional unknown bases (N) into the assembly because the length of each stretch was estimated based on the tobacco genome.

Repeat content

More Tnt1 retrotransposons are found in N. tomentosiformis than in N. sylvestris, (7.39% and 3.98% respectively), which apparently contradicts previous reports [42–44]. This finding could be caused by the mislabeling of novel N. tomentosiformis repetitive elements obtained by RepeatScout as Tnt1. The amounts of Tnt2 and Tto1 repetitive elements are higher in N. sylvestris than in N. tomentosiformis and this finding agrees with previous studies. Furthermore, as reported previously [45], we also observed a higher proportion of NicCL3 (0.75% versus 0.14%) and NicCL7/30 (5.57% versus 2.94%) repetitive DNA elements in N. tomentosiformis than in N. sylvestris.

Genetic markers

The 2,363 tobacco SSR markers reported previously [13] were mapped to both genome assemblies. The number of uniquely mapped markers on each genome was then compared with the results of the PCR amplification tests performed in N. sylvestris and N. tomentosiformis, in order to assign an origin to them when creating the tobacco genetic map (Additional file 4). Sixty-five percent of the SSR markers that amplified only in N. sylvestris mapped only to the N. sylvestris genome; 7% mapped to both genomes. Similarly, 65% of the SSR markers that amplified only in N. tomentosiformis mapped only to N. tomentosiformis; 15% mapped to both N. sylvestris and N. tomentosiformis. About a third of the tobacco SSR markers could not be mapped. This can be expected, because the current draft genome assemblies are likely to fail assembling in regions with simple repeats such as the ones found in SSR markers. If this is the case, a primer pair will match to two different sequences.

Of the 173 SSR markers present in the N. acuminata genetic map (Additional file 5), 128 (74%) of them could be mapped to the N. sylvestris genome assembly. This number is the sum of the 75 SSRs of the N. acuminata map found in the N. sylvestris assembly, the 50 SSRs of the N. acuminata map found in the N. sylvestris and N. tomentosiformis assemblies, the single SSR of the N. acuminata and N. tomentosiformis maps found in the N. sylvestris assembly, and the 2 SSRs of the N. acuminata and N. tomentosiformis maps found in the N. sylvestris and N. tomentosiformis assemblies (Additional file 6). Similarly, of the 221 SSR markers present in the N. tomentosiformis genetic map (Additional file 7), 173 (78%) could be mapped to the N. tomentosiformis genome assembly (Additional file 6). In addition, 706 SSR markers not present on the existing genetic maps could be mapped to the N. sylvestris genome assembly, 605 mapped to the N. tomentosiformis genome assembly, and 174 mapped to both.

Of the 134 COSII markers present in the N. acuminata genetic map, 45 (34%) could be mapped to the N. sylvestris genome assembly (Additional file 8). Similarly, of the 262 COSII markers in the N. tomentosiformis genetic map, 81 (31%) could be mapped to the N. tomentosiformis genome assembly (Additional file 8). Using the same method, 736 of the 879 COSII markers (84%) on the expen2000 tomato genetic map could be found; 718 of them mapped to the expected chromosome. In addition, 68 COSII markers not present on the existing genetic maps could be mapped to the N. sylvestris genome assembly, 78 mapped to the N. tomentosiformis genome assembly, and 226 mapped to both.

The low numbers of COSII markers that could be mapped to the N. sylvestris and N. tomentosiformis assemblies, despite the good results that were obtained using the same method on the tomato map, could be due to the current fragmented state of the assemblies, or because the COSII marker primers are not adapted for Nicotiana species.

Transcriptome assembly

Transcript and protein quality

The assembled reference transcriptome was assessed for completeness and accuracy by mapping the transcripts to the UniProt reference plant sequence databases. The number of sequences for both the transcripts and the unique genes from which the transcripts are derived that could be mapped was similar for N. sylvestris and N. tomentosiformis (Figure 1). For N. sylvestris and N. tomentosiformis, 58.6% and 60.5% of transcripts, respectively, had significant ORFs with a length equal to or longer than 100 amino acids. The majority, 82.2% for N. sylvestris and 81.9% for N. tomentosiformis, had a homologous sequence in the UniProt Knowledgebase. Approximately a third of these peptide sequences, 37.2% in N. sylvestris and 36.5% in N. tomentosiformis, had hits in Swiss-Prot, the annotated subset of UniProt. The BLAST alignments show that while the coverage of the predicted ORFs by the reference sequences is generally high (Figure 2) and comparable between the species, the coverage of the reference sequence by the predicted ORFs is often partial, indicating that these ORFs are likely to be incomplete.

Functional comparison to other species

We used the OrthoMCL software [54, 55] to define clusters of orthologous and paralogous genes between N. sylvestris and N. tomentosiformis, as well as tomato, another representative of the Solanaceae family, and Arabidopsis as a representative of the eudicots (Figure 3). While a large number of sequences are shared between all the species (7,098), many are specific to Solanaceae (2,790). A very high number of sequences are only observed in the Nicotiana species (3,614), with several hundred gene clusters being specific to N. sylvestris and N. tomentosiformis. These sequences may be artifacts that are the result of incomplete transcripts not clustering correctly, rather than actual novel protein families that evolved since the split of the species.

At the tissue level (Figure 4), the vast majority of gene clusters are shared. As far as the number of clusters is concerned, flowers had the most diverse transcriptome; flowers also contain a large number of transcripts (3,510 for N. sylvestris and 3,387 for N. tomentosiformis) not found in root or leaf tissues. The number of tissue-specific clusters is very low (<20 for any of the tissues in either of the species); this number reflects the 'noise level' of the merging process because in choosing representative transcripts while merging of the tissue transcriptomes, a different set of exons may have been chosen, and the tissue sequences may not match the representative in the merged transcriptome.

Functional annotation

Over 7,000 proteins could be annotated with a three-digit EC number using the EFICAz tool (7,554 in N. sylvestris, 7,397 in N. tomentosiformis), of which over 4,000 were assigned with high confidence (4,061 in N. sylvestris and 4,008 in N. tomentosiformis). This implies that just less than 20% of the predicted proteome of the two species (19.4% and 19.1% for N. sylvestris and N. tomentosiformis) has enzymatic function. Just over 4,000 and over 3,000 four-digit EC numbers could be assigned to predicted proteins. Although the number of unique four-digit EC numbers is comparatively small (635 for both N. sylvestris and N. tomentosiformis), this information can still be used to generate molecular pathway databases.

Approximately half of all the proteins were annotated with at least one GO term by the InterProScan software (57.8% for N. sylvestris and 57.4% for N. tomentosiformis); close to 50,000 biological process tags were assigned and slightly more than 20,000 molecular functions were assigned to just under 20,000 unique proteins. GO term enrichment was analyzed using the GOStats package [56]. Enrichment was tested by comparing the GO term complement for each species against the background of the pooled set of GO terms from both organisms. We see only small and not highly significant changes in gene composition. For N. sylvestris, the defense response function is overrepresented; in N. tomentosiformis we observe an enrichment of core metabolic functions as well as protein phosphorylation. The phenotypic differences between the species are thus likely to be regulatory rather than due to the loss or gains of new genes. A caveat exists in that the absence of a gene in the genome of one species does not guarantee that the gene does not exist: it is possible that the section of the genome containing the gene simply has not been covered by our current efforts and that further sequencing will identify these. The results of GO term enrichment analysis are shown in Additional file 11.

Heavy-metal transport

A non-exhaustive list of gene copies that may be involved in cadmium/zinc (Cd/Zn) accumulation (Figure 5) in Nicotiana leaves is shown in Additional file 12. The corresponding transcripts in root, leaf and flower are depicted. The expression data resulting from the hybridization of specific Affymetrix probes (100% matches with the targeted sequences) with leaf RNA isolated from N. sylvestris and N. tomentosiformis provided data similar to fragments per kilobase of transcript per million mapped reads (FPKM) expression data. The results show that the design of the Affymetrix exon probes is suitable for the analyses of gene expression in both tobacco ancestors, N. sylvestris and N. tomentosiformis.

Based on sequence and expression data analogies with corresponding Arabidopsis thaliana gene data, two Nicotiana iron transport-related sequences belonging to the IRT family were identified and named NsylIRT1, NtomIRT1 and NsylIRT2, NtomIRT2 corresponding to Arabidopsis IRT1 and IRT2. Both of the A. thaliana genes are expressed in the roots and are involved in Zn/Cd uptake [57], although IRT1 is more selective for iron [58]. Interestingly, IRT1 and IRT2 are expressed in N. sylvestris roots but not in N. tomentosiformis roots, suggesting that one or more other genes, possibly belonging to the ZIP family, function for Zn and iron uptake in N. tomentosiformis [59]. Conversely, the potential Nicotiana orthologs (NtomIRT3 and NsylIRT3) of AtIRT3 are not expressed in the roots (Additional file 12), although AtIRT3 is expressed in Arabidopsis roots, where it is involved in Zn and iron transport [60]. Interestingly, NsylIRT3 and NtomIRT3 transcripts are more abundant in flower tissues most likely for the redistribution of Zn and Fe. The function of Nicotiana IRT3 is possibly closer to the Zrt/IRT-like protein AtZIP4, which is highly expressed in anther and pollen [61], where it is suspected to play a role in Zn redistribution in flowers ([62] and references therein). Thus, Zn and iron uptake is likely driven by AtIRT1 and AtIRT2 orthologous proteins in N. sylvestris, whereas another gene(s) is likely to perform this function in N. tomentosiformis.

The P_1B-type ATPases, known as heavy metal ATPases (HMAs), play important roles in metal transport in plants. In Arabidopsis, AtHMA3 (OsHMA3 in rice) is localized in the tonoplast membrane, where it plays a major role in detoxifying Zn and Cd through vacuolar sequestration. AtHMA3 is recognized as the major locus responsible for the variation in leaf Cd accumulation of A. thaliana accessions. AtHMA2 (OsHMA2 in rice) and AtHMA4 (likely OsHMA9 in rice), are localized in the plasma membrane and are expressed in the tissues that surround the vascular vessels of roots, where they function in Zn and Cd efflux from cells [63]. In N. sylvestris, N. tomentosiformis and Solanum lycopersicum (Solyc07g009130) genomes, only one HMA gene orthologous to the sub-cluster formed by AtHMA2, AtHMA3 and AtHMA4 in A. thaliana is present. This suggests a strong evolutionary divergence between Brassicaceae-Poaceae and Solanaceae. The FPKM expression data show major expression of Nicotiana HMA in the root tissues, suggesting that it has functions that are similar to those of AtHMA2, AtHMA3 and AtHMA4, and is more involved in Zn/Co/Cd/Pb translocation from root to shoot than in vacuolar sequestration.

The long-distance root-to-shoot transport of Cd/Zn can be driven by phytochelatins or nicotianamine. Therefore, the key genes that may affect Cd/Zn accumulation in leaves are phytochelatin synthases (PCS) and nicotianamine synthetases (NS) [64, 65]. The orthologous genes (PCS, NS1 and NS2) identified in N. sylvestris and N. tomentosiformis exhibit similar expression profiles in the root, leaf and flower tissues, suggesting that transport in vascular tissues is similar in both Nicotiana species.

Genes orthologous to the ABC transporters that are involved in Cd transport in A. thaliana, such as AtPDR8 (PDR) [66] and AtATM3 (ATM) [67], are found in both the N. sylvestris and N. tomentosiformis genome. Their expression profiles are similar in both Nicotiana species and close to their expression profiles in Arabidopsis, suggesting that these genes have similar functions in root, leaf and flower in both species. ABC proteins related to the multidrug resistance-associated protein (MRP) family have been already described to be involved in Cd transport and sequestration [68, 69] Although the precise cellular function of one of the MRP family members in N. tabacum, NtMRP4, has not yet been determined, silencing NtMRP4 [70] resulted in Cd reduction in leaves under field conditions. The corresponding gene is expressed in both N. sylvestris and N. tomentosiformis, suggesting that it has similar functions in both plants.

Other genes that play a role in Cd accumulation into vacuoles belong to the cation proton exchanger (CAX) family. Overexpression of AtCAX2 and AtCAX4 in tobacco resulted in Cd accumulation in the roots and a subsequent decrease in the shoots [71]. Four genes that clustered with AtCAX2 and not AtCAX4 (Additional file 13) were identified in N. sylvestris and N. tomentosiformis, suggesting that tobacco CAX gene products orthologous to AtCAX2 and not AtCAX4 may play roles in Cd sequestration in Nicotiana species. The expression profiles of the four genes are similar in both N. sylvestris and N. tomentosiformis, indicating that these genes play identical functions in both plants.

Alkaloid metabolism

The key genes involved in the synthesis of nicotine and nornicotine alkaloids in Nicotiana leaves (Figure 6) are listed in Additional file 14 and the corresponding transcripts in root, leaf and flower are shown. The expression data obtained from the hybridization of specific Affymetrix probes (100% match with the targeted sequences) with leaf RNA isolated from N. sylvestris and N. tomentosiformis provided data similar to FPKM expression, except for four N. tomentosiformis genes - NtomQPT1, NtomBBL3, NtomNND1 and NtomNND2. However, these four genes were found to be expressed in the leaf of N. tomentosiformis plants subjected to RNA-seq analyses. The plants that were used for the RNA-seq analyses were fully mature compared with the young plantlets that were used for the Tobacco Exon Array hybridization, which may indicate that the four genes are more highly expressed in mature leaves than in the primary leaves, suggesting that these genes may possibly affect the alkaloid pathway. Similar to the Cd genes described above, this type of comparison confirms that the design of the Affymetrix exon probes is suitable for the analyses of gene expression in both N. sylvestris and N. tomentosiformis.

The higher accumulation of nicotine in N. sylvestris compared with N. tomentosiformis is because of the relatively large deletion that encompasses the NIC2 locus of N. tomentosiformis. Therefore, the low-nicotine phenotype is often associated with nic2 mutations [72]. In nic1nic2 mutant roots, BBL transcripts are strongly reduced [73], attesting that berberine bridge enzyme-like (BBL) genes are regulated by the NIC loci in the roots. Our data confirm that BBL1 and BBL3 are particularly expressed in the roots of both Nicotiana species. However, no large differences in transcript levels were found, possibly suggesting that BBL gene regulation is not as different as suspected between N. sylvestris and N. tomentosiformis, and the effect of the nic2 deletion is apparent somewhere else within the nicotine biosynthesis pathway. In this context, our data show that the expression of a large set of genes involved in nicotine biosynthesis, for example, L-aspartate oxidase (AO), quinolinate synthase (QS), quinolinate phosphoribosyltransferase (QPT), and putrecine N-methyltransferase (PMT), are strongly up-regulated in the roots of N. sylvestris compared with N. tomentosiformis; indeed, PMT expression is not detected in the roots of N. tomentosiformis. Four different PMT genes have been found in N. tabacum [74] and, based on sequence analogy, three of them [75] likely originate from N. sylvestris (Additional file 15). Surprisingly, the two copies of PMT that are present in N. tomentosiformis are similar to only one PMT gene in N. tabacum (NtPMT-2). This finding suggests that because of the lack of the three other PMT copies in N. tomentosiformis, the full pathway for nicotine synthesis is certainly different in N. tomentosiformis than in N. sylvestris, which has three PMT copies that are related to N. tabacum, NtPMT-1, -3 and -4 (Additional file 15).

The up-regulation of PMTs, AO and QS in N. sylvestris compared with N. tomentosiformis attests that the early steps in the pathway that lead to the synthesis of nicotinic acid are also particularly active in N. sylvestris and certainly play a major role in nicotine synthesis. Recent data reported by Shoji and Hashimoto [76] suggest that tobacco MYC2 regulates PMT-2 and QPT-2 by interacting with specific promoter regions. It is therefore tempting to speculate that regulation occurs differently via MYC2 in N. sylvestris and N. tomentosiformis. Conversely, because AO and QS are located in the plastids and are involved in NAD synthesis from aspartate via quinolinic acid [77], they are likely regulated via nuclear cross-talk that is possibly more active in N. sylvestris than in N. tomentosiformis (Additional file 14).

In species of the Nicotiana genus, the conversion of nicotine to nornicotine, which is the precursor of the tobacco nitrosamine N'-nitrosonornicotine, is mediated by nicotine N-demethylase enzymes (NND) encoded by the CYP82E subfamily of cytochrome P450 genes. Four genes from this gene family are reported to be distributed in the N. sylvestris and N. tomentosiformis genomes. CYP82E4 is the dominant factor in senescence-inducible nornicotine production, whereas CYP82E5v2 is involved in nicotine conversion in the green leaves; both of them are found in N. tomentosiformis, along with CYP82E3. In N. sylvestris, one such gene, CYP82E2, has been found [78]. Searches in both these Nicotiana genomes revealed that N. sylvestris has five related genes, and N. tomentosiformis has four. A comparison of the phylogenetic trees (Additional file 16) confirms that three of the N. tomentosiformis genes are related to the N. tabacum CYP82E3, CYP82E4 and CYP82E5 genes and that one of the N. sylvestris genes corresponds to N. tabacum CYP82E10 [79]. The data presented in Additional file 14 and Additional file 16 show that NtomNND-1 is evolutionarily close to one copy of CYP82E4and highly expressed in flowers, whereas its expression in leaves is not supported by Affymetrix data. To our knowledge, the high expression of a nicotine demethylase gene in flowers has not yet been described; the gene product possibly plays a role in protection against insects. Conversely, the NsylNND-1 that is evolutionary close to the N. tabacum CYP82E10 is highly expressed in roots, confirming the findings of an earlier study [79]. The high expression of the three N. tomentosiformis genes related to the N. tabacum CYP82E3, CYP82E4 and CYP82E5 genes suggests that N. tomentosiformis is globally a more active producer of nornicotine than N. sylvestris, which is the opposite of what was found for nicotine synthesis (see above).

Illumina sequencing

Young leaves, roots and flowers of N. sylvestris (USNGC TW136, PI555569) and N. tomentosiformis (USNGC TW142, PI555572) grown in a greenhouse were collected. DNA extraction was performed using Qiagen DNAeasy Plant Maxi Kit (Qiagen, Düsseldorf, Germany) from fresh leaves. RNA extraction was performed using the Qiagen RNAeasy Mini Kit (Qiagen).

Short insert 'paired-end' libraries were prepared using the Illumina TruSeq DNA Sample Preparation Kit version 2 according to the manufacturer's instructions, or with few modifications if prepared by Fasteris. For Fasteris, 2.1 mg of genomic DNA was broken using BioRuptor (Diagenode, Liège, Belgium); ends were repaired using Klenow and polynucleotide kinase, and then Fasteris-modified adapters were ligated to the inserts. After size selection on agarose gel, the libraries were amplified by ten PCR cycles, and then purified and quantified.

Long insert 'mate-pair' libraries were prepared using the Illumina Mate Pair Library Prep Kit version 2 according to the manufacturer's instructions, or using a Fasteris-developed protocol in which 10 mg of genomic DNA were broken into fragments of approximately 2 to 5 kb using Covaris (KBioSciences, Herts, UK) and purified on 0.7% agarose gel to recover fragments of 3 kb and 5 kb. After end repair, a Fasteris-designed spacer was ligated and the fragments were circularized. Non-circular fragments were eliminated and then the DNA was broken using Covaris to generate fragments of 400 bp, which were end repaired, ligated with Illumina adapters, purified on agarose gel and amplified by PCR for 12 cycles.

RNA-seq libraries were constructed using Illumina's TruSeq RNA Sample prep Kit protocol according to the manufacturer's instructions. All the libraries (short insert 'paired-end', long insert 'mate-pair', and RNA-seq) were sequenced on an Illumina HiSeq-2000 using version 3 chemistry and flow-cells with runs of 2 × 100 bases. Base calling and sample demultiplexing were performed using Illumina's HiSeq Control Software and the CASAVA pipeline. The data for the N. sylvestris and N. tomentosiformis RNA-seq triplicates have been uploaded to the EBI Sequence Read Archive under accession numbers ERP002501 and ERP002502, respectively.

Genome size estimation

We estimated the genome size of N. sylvestris and N. tomentosiformis using the 31-mer depth distribution of all the non-overlapping paired-end libraries, as described previously [36, 50–52]. Briefly, the genome size is obtained by dividing the total number of 31-mers considered to be error-free by their most frequent depth of coverage.

Genome assembly

The raw DNA reads from N. sylvestris and N. tomentosiformis were preprocessed by first trimming 3' bases with qualities lower than 30, and then discarding reads shorter than 50 bases or with less than 90% of the bases with qualities lower than 30. The paired-end libraries with insert sizes shorter than 200 bases were further preprocessed using FLASH [80] to merge the paired-end reads into extended single reads.

The paired and single reads from the paired-end libraries were then assembled into contigs using SOAPdenovo [81] with a k-mer of 63, and the paired reads from paired-end and mate-pair libraries were used for scaffolding by increasing library size. To improve scaffolding, mate-pair libraries from closely related Nicotiana species were also used. Gaps that resulted from the scaffolding were closed using GapCloser and all sequences shorter than 200 bases were discarded from the final assemblies.

Superscaffolding using the tobacco WGP™ physical map was possible because it is based on sequencing tags, and the origin of the WGP contigs have been annotated. Briefly, WGP tags of S or T origin were mapped to the N. sylvestris or N. tomentosiformis sequences, respectively. Superscaffolds were created when two or more sequences could be anchored and oriented unambiguously to a WGP contig. The N. sylvestris and N. tomentosiformis genome assemblies have been submitted to GenBank BioProjects PRJNA182500 and PRJNA182501, respectively. The N. sylvestris whole genome shotgun project has been deposited at DDBJ/EMBL/GenBank under the accession ASAF00000000. The version described in this paper is version ASAF01000000. The N. tomentosiformis whole genome shotgun project has been deposited at DDBJ/EMBL/GenBank under the accession ASAG00000000. The version described in this paper is version ASAG01000000. The raw sequencing data used for the assemblies of N. sylvestris and N. tomentosiformis genomes have been submitted to the EBI Sequence Read Archive under accession numbers ERP002501 and ERP002502.

Repeat content estimation

The repeat content of the N. sylvestris and N. tomentosiformis genome assemblies were estimated using RepeatMasker [82] with the eudicot repeat library available from the Sol Genomics Network, the TIGR Solanaceae repeat library, and RepeatScout [83] libraries created using sequences of at least 200 kb from the draft genome assemblies of N. sylvestris and N. tomentosiformis. Classification of the repeat types was done using the NCBI BLASTN [84–86] hits to known repeat elements.

Genetic markers

PCR primers for the SSR markers have been reported previously [13] and the COSII makers from Sol Genomics Network were mapped to the draft assembly genomes of N. sylvestris and N. tomentosiformis using LAST [87]. Only the primer pairs that could be mapped with at least 95% identity and that yielded a unique PCR product were retained.

Pathway gene identification and quantification

Genomic regions containing genes that potentially encode proteins from the selected pathways were identified by mapping homologous proteins from other species to the genome assemblies using BLAT [88] and manually curating the hits. Probes from the Tobacco Exon Array [49] were selected by mapping them to the identified genome regions using LAST [87] and retaining only perfect matches that could be mapped uniquely. Quantification of gene expression was obtained by summing the Cufflinks [89] FPKM values of the transcripts that overlapped the identified genome regions.

De novo transcriptome assembly

All the reads were preprocessed to clip the overrepresented sequences reported by FastQC [90]. After clipping, the 3' ends of the reads were quality trimmed with a quality threshold of 20 and artifacts were removed. Finally, reads of at least 50 nucleotides with at least 75% nucleotides of quality 20 or more were kept. The clipping, trimming and filtering were performed using the fastx toolkit [91]. Transcripts were assembled using the Trinity de novo assembly pipeline [92]; the peptide prediction program contained within this software suite was used to predict peptides from the assembled transcripts.

Transcriptome assembly was performed using the 'Tuxedo' suite of tools. Reads were mapped to the appropriate genome assembly using the Bowtie2/Tophat2 [93, 94] pipeline with the default parameters. Transcript generation was performed using the Cufflinks tools [89] and merged using Cuffmerge. A representative set of transcript sequences was generated using the 'gtf_to_fasta' component of Cufflinks.

Transcript and protein quality

The ORF finding utility included in the Trinity software package (version January 2012) was used to find ORFs in the inferred transcripts. Candidate peptide sequences were culled at a minimum length of 100 amino acids. The search for sequences homologous to the ORFs was performed using BLAST [84], with the UniProt Knowledgebase and the Swiss-Prot subset as reference databases. A reasonably stringent e-value cutoff of 1E-30 was used and only one hit was retained for each sequence. To determine which of the best hits were mutual, a reverse search (reference database against transcripts) was also performed using the same parameters.

Functional comparison to other species

Orthologous and paralogous genes between our sequences and those from other species were clustered using OrthoMCL [54, 55]. To ensure comparability, we used the same ORF finding software on the Arabidopsis (TAIR10 sequence) and tomato sequences (ITAG version 2.3) to derive peptide sequences and then used only sequences of 100 amino acids or longer. An all-against-all sequence search was performed using BLAST [84] (version 2.2.23+) with default parameters and the results of this search were used as the input to OrthoMCL, which was run using the default parameter set (percentMatchCutoff = 50; e-valueExponentCutoff = -5). The OrthoMCL protein group output files were further processed using in-house Python scripts, and visualized in R [95] as a Venn diagram using the CRAN package 'VennDiagram' [96].

Functional annotation

The EFICAz^2.5 software [97–99] was used to predict EC numbers for the protein sequences predicted from the transcripts of the pooled tissue samples. The InterProScan software [100], version 4.8 was used to assign GO terms to the protein sequences.