Transcriptome-enabled marker discovery and mapping of plastochron-related genes in Petunia spp.

Petunia (Petunia × hybrida), derived from a hybrid between P. axillaris and P. integrifolia, is one of the most economically important bedding plant crops and Petunia spp. serve as model systems for investigating the mechanisms underlying diverse mating systems and pollination syndromes. In addition, we have previously described genetic variation and quantitative trait loci (QTL) related to petunia development rate and morphology, which represent important breeding targets for the floriculture industry to improve crop production and performance. Despite the importance of petunia as a crop, the floriculture industry has been slow to adopt marker assisted selection to facilitate breeding strategies and there remains a limited availability of sequences and molecular markers from the genus compared to other economically important members of the Solanaceae family such as tomato, potato and pepper. Here we report the de novo assembly, annotation and characterization of transcriptomes from P. axillaris, P. exserta and P. integrifolia. Each transcriptome assembly was derived from five tissue libraries (callus, 3-week old seedlings, shoot apices, flowers of mixed developmental stages, and trichomes). A total of 74,573, 54,913, and 104,739 assembled transcripts were recovered from P. axillaris, P. exserta and P. integrifolia, respectively and following removal of multiple isoforms, 32,994 P. axillaris, 30,225 P. exserta, and 33,540 P. integrifolia high quality representative transcripts were extracted for annotation and expression analysis. The transcriptome data was mined for single nucleotide polymorphisms (SNP) and simple sequence repeat (SSR) markers, yielding 89,007 high quality SNPs and 2949 SSRs, respectively. 15,701 SNPs were computationally converted into user-friendly cleaved amplified polymorphic sequence (CAPS) markers and a subset of SNP and CAPS markers were experimentally verified. CAPS markers developed from plastochron-related homologous transcripts from P. axillaris were mapped in an interspecific Petunia population and evaluated for co-localization with QTL for development rate. The high quality of the three Petunia spp. transcriptomes coupled with the utility of the SNP data will serve as a resource for further exploration of genetic diversity within the genus and will facilitate efforts to develop genetic and physical maps to aid the identification of QTL associated with traits of interest.


Background
The genus Petunia resides within the Solanaceae family and contains 20 species and subspecies that are native to South America [1]. Petunia × hybrida (petunia) is an important ornamental crop plant and represents a hybrid species derived in the nineteenth century from a cross between P. axillaris and P. integrifolia [2]. Subsequent breeding has introgressed traits from additional Petunia spp. that, together with natural variation resulting from mutations in key genes, have contributed to the wide diversity in plant and floral morphology and flower color that exists within the pool of commercially available germplasm [2][3][4][5][6][7]. In cool climates in the Northern Hemisphere, petunia is often produced in greenhouses during the winter months for distribution to spring markets once it reaches an optimal size and begins to flower [8,9]. Therefore, a high percentage of the cost of crop production is related to energy consumption and growers are often faced with the dilemma of either reducing greenhouse temperatures, thereby extending the growing time of the crop and incurring increased labor costs, or elevating the growing temperature and increasing energy costs but reducing the duration of crop growth [8,9]. Thus, understanding the factors that impact crop timing traits may facilitate the selection of petunia varieties with an increased rate of vegetative node formation (development rate) at either optimal or suboptimal growing temperatures, or varieties that initiate flowering following emergence of fewer leaf nodes. We have previously documented that accessions of P. axillaris and P. integrifolia possess increased development rate when compared to a diverse pool of commercial petunia germplasm, suggesting genetic variation for this trait within the genus [10]. This was confirmed in an interspecific F 2 population of a cross between P. axillaris and P. integrifolia that identified three quantitative trait loci (QTL) on chromosomes 1, 2 and 5 that affected development rate and explained 34 % of the observed variation [11]. The molecular basis underlying these QTL remains to be identified.
The genetic determinants of development rate, often referred to as plastochron, are multifaceted, complex and not fully understood but are, at least in part, linked to hormonal control of meristem size and activity. For example, transgenic tobacco (Nicotiana tabacum) plants with increased cytokinin oxidase activity and a concomitant reduction in cytokinin levels displayed reduced meristem size and delayed plastochron when compared to wild type [12]. Similarly, characterization of Arabidopsis mutants with reduced auxin levels and disrupted auxin transport also influence plastochron [13][14][15]. In addition, mutations at the plastochron 1 and plastochron 2 loci of rice, which encode a cytochrome P450 of unknown function and a MEI2-like RNA binding protein homolog, respectively, also influence development rate but do so independently of each other [16,17]. In Arabidopsis, the SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) transcription factors SPL9 and SPL15 act redundantly to influence plastochron and over-expression of miR156, which targets multiple SPLs, shortens plastochron [18]. Analogous to the relationship of the plastochron 1 and 2 loci of rice [17], the SPL/miR156 regulatory module acts independently of CYP78A5/KLUH, which is a putative ortholog of rice plastochron 1 [18]. The involvement of the miRNA pathway in influencing plastochron is further supported by the characterization of the serrate and altered meristem program (amp1) mutants of Arabidopsis, which display reduced and increased rates of leaf initiation, respectively [19][20][21][22]. SERRATE encodes a zinc finger protein required for miRNA biogenesis and RNA splicing while AMP1 associates with ARGONAUTE1 at the endoplasmic reticulum and is required for translation inhibition through the exclusion of miRNA target mRNAs from polysomes [19,20,23]. Furthermore, mutations in AMP1 homologs in maize and rice confer similar pleiotropic phenotypes to those exhibited by Arabidopsis amp1 mutants, including altered plastochron [24,25]. Together, these data suggest complex regulation of plastochron that involves different regulatory modules, including hormone and miRNA pathways.
The development of next generation sequencing technology has revolutionized biology and in particular, transcriptome sequencing provides a cost effective strategy for generating sequence and expression information from the gene space of non-model organisms or from species with large complex genomes [26,27]. In plants, transcriptome sequencing has facilitated gene discovery, the development of molecular markers and large scale analyses of genetic variation [28][29][30][31][32][33]. Despite the economic and biological importance of petunia, genomic information and molecular marker resources for this genus are limited [11,[34][35][36][37][38], single nucleotide polymorphism (SNP) markers are currently unavailable and marker assisted selection is rarely utilized. In addition, although transcriptome resources are available for petunia, they are not extensive and most often are derived from the cultivated species P. hybrida or from specialized tissue types [39][40][41]. Herein, de novo assembly of reference transcriptomes of P. axillaris, P. integrifolia and P. exserta, derived from paired-end RNAseq analysis of five diverse tissues (callus, seedling, shoot apices, flowers and trichomes) is reported. Tissue types were selected to attempt to maximize the number of transcripts recovered while generating resources for studying traits of interest related to development and metabolism. These resources were utilized to develop a set of SNP, cleaved amplified polymorphic sequence (CAPS) and simple sequence repeat (SSR) markers that will facilitate QTL mapping and gene discovery for multiple traits within the genus, including those associated with development rate.

Results and discussion
Transcriptome assembly and annotation Transcriptome sequencing of five tissue libraries, including callus, flowers, shoot apex, seedlings, and trichomes, from P. axillaris, P. exserta, and P. integrifolia yielded between~248 and 294 M 100 nt reads, of which greater than 94 % passed quality and trimming filters (Table 1). A two-step de novo assembly strategy (see Additional file 1: Figure S1) modified from [42] resulted in 74,573, 54,913, and 104,739 transcripts ≥ 500 bp respectively, for P. axillaris, P. exserta, and P. integrifolia. The two-step assembly strategy was employed to eliminate redundant reads in the flower and trichome libraries already present in the callus, shoot apex, and seedling libraries to aid the quality of the assembly.
The N50 value for each assembly was greater than 1950 and average transcript sizes of 1714 for P. axillaris, 1624 for P. exserta, and 1646 for P. integrifolia were obtained. The distribution of transcript sizes follows the same trend in each of the three species with the largest number of transcripts falling within the size bins of 500-1000 bp, and between 1001 and 1500 bp (Fig. 1a). CEGMA analysis [43] revealed that full length copies of >91 % of the highly conserved eukaryotic genes are present in each of the transcriptome assemblies while partial sequences are present for almost 100 % of these genes (Fig. 1b). These data are similar to those reported for other plant transcriptome assemblies [44,45]. All high quality raw reads from each of the five libraries were also mapped back to their respective assembly with Bowtie and Tophat [46] generating mapping rates of 87.2 % for P. axillaris, 86.8 % for P. exserta, and 78.3 % for P. integrifolia.
For P. axillaris and P. exserta, the majority of representative transcripts comprised a single isoform whereas approximately 60 % of P. integrifolia transcripts possessed multiple isoforms (Table 2; Additional file 1: Figure S2). Where multiple isoforms were identified, the median number was two for P. axillaris and P. exserta, and three for P. integrifolia. The longest isoform was selected as the representative transcript yielding a total of 32,994, 30,225 and 33,540 representative transcripts in P. axillaris, P. exserta, and P. integrifolia, respectively. These transcripts covered between~47.5 and 53.1 Mbp of the transcriptome space ( Table 2). The increased number of transcripts retrieved from the P. integrifolia assembly together with a higher number of transcripts with multiple isoforms, is likely the result of widespread heterozygosity within this species due to self-incompatibility [47].
Representative transcripts from each assembly were annotated using BLASTX searches against five publically available databases (Additional file 2: Dataset S1). Adopting quality thresholds of ≥30 % coverage and ≥70 % identity resulted in an annotation rate of approximately 70 % (Fig. 1c), with the largest number of annotations retrieved from the RefSeq and NCBI nonredundant databases. Similar rates of annotation were reported for transcriptome assemblies of chickpea and red clover [48,49]. In addition, of the 1944 predicted Petunia proteins available in GenBank,~89 % are present in the P. axillaris, P. exserta, and P. integrifolia transcriptome assemblies.
Open reading frames (ORFs) were extracted from the representative transcripts and the predicted protein sequences were searched for orthologous gene clusters using OrthoMCL [50]. Among all comparisons, approximately 21,000 orthologous clusters were identified, with P. axillaris, P. exserta, and P. integrifolia sharing 13,747 clusters ( Fig. 2; Additional file 3: Dataset S2). A larger number of orthologous clusters were found between P. axillaris and P. exserta, likely due to their closer phylogenetic relationship [1]. While P. integrifolia was found to house over double the number of unique clusters ( Fig. 2; Supplemental Dataset 2). In addition, Gene Ontology (GO) annotations (Additional file 4: Dataset S3) of the representative transcripts showed highly equivalent representation of biological process categories The number and percentage of total raw reads included in each transcriptome assembly after quality filters were applied Fig. 1 Quality metrics of the Petunia spp. transcriptome assemblies. a Size distribution of assembled transcripts. b CEGMA completeness assessment of the transcriptome assemblies. c Percentage of P. axillaris, P. exserta, and P. integrifolia unigenes with assigned functional annotations from UniRef100, TAIR10, RefSeq, the Pfam domain database, and NCBI GenBank non-redundant protein set. In addition, the total percentage of annotated unigenes per species is presented among the three species, indicating the uniformity of the transcriptome assemblies and their subsequent annotation (Fig. 3).

SNP detection and characterization
In the absence of whole genome sequences, comparative transcriptome analysis has proven utility for developing SNP markers [30,31]. To facilitate genetic analysis within the genus, the three Petunia spp. transcriptomes were utilized for SNP discovery using the Genome Analysis Toolkit (GATK; see Additional file 1: Figure S3) [51]. When utilizing RNA-seq data for SNP discovery, removal of duplicate reads increases SNP detection sensitivity and specificity [52], thus this strategy was adopted resulting in utilization of 62.7 % of P. axillaris, 51.3 % of P. exserta and 60.8 % of P. integrifolia reads that uniquely mapped to the P. axillaris transcriptome (data not shown). Depth of sequence coverage increases the reliability of SNP detection [52] and utilization of a minimum threshold for read coverage depth of 10 resulted in mapping of 93.7 %, 81.2 % and 73.3 % of the reads for P. axillaris, P. exserta, and P. integrifolia, respectively (Table 3).
Overall, there were 814,408 SNPs between the three species. Of these, 105,645 were located in 5ʹ untranslated regions (UTRs), 481,289 were located within the coding sequence (CDS), and 217,178 were located in 3ʹUTRs. SNP frequency was calculated by dividing the total length of the reference transcriptome by the total number of SNPs. When only considering the length of 5ʹUTRs, CDS, and 3ʹUTRs in transcripts (not the length of the genomic regions these transcripts spanned), the SNP frequencies were 1 SNP/69 bp, 1 SNP/61 bp and 1 SNP/62 bp in these regions, respectively. The overall SNP frequency was 1 SNP/63 bp.
After filtering, among the 32,994 representative P. axillaris transcripts (unigenes), we identified SNPs between P. axillaris and either P. exserta or P. integrifolia in 20,606 unigenes. Gene Ontology (GO) annotation revealed that among all unigenes, 22,535 (68.3 %) contained GO terms, while 16,787 (80.8 %) of the SNPcontaining unigenes were assigned with one or more GO ID (Additional file 1: Figure S4). In general, the distribution of GO terms was very similar between all unigenes and those containing SNPs. KEGG Pathway analysis was carried out to determine functional categorization of unigenes containing SNPs. A total of   Table 4). These unigenes were assigned to 322 KEGG pathways. Overall, the SNP frequency for genes involved in metabolism and organismal systems super pathways was lower than those involved in genetic information processing and environmental information processing.
In total, 89,007 SNP positions were detected among P. axillaris, P. exserta and P. integrifolia. Of these, 8868 (10.0 %) were polymorphic across all three species combinations. This calculation was based on positions for which at least one species, P. integrifolia or P. exserta, had sufficient depth of coverage (10X). After removing low depth of coverage SNPs, 73,193 SNPs remained between P. axillaris and P. integrifolia, 25,847 SNPs remained between P. axillaris and P. exserta, and 79,438 SNPs remained between P. exserta and P. integrifolia. (Fig. 4a).
SNPs were further classified based on their location; CDS or UTR and their zygosity. Altogether, only 10 homozygous SNPs were polymorphic across all three species combinations (Fig. 4b). Between P. axillaris and P. exserta, 24,528 SNPs were homozygous, which comprised 94.9 % of all SNPs. There were 13,512 homozygous SNPs, 17 % of the total SNP positions, between P. exserta and P. integrifolia. Between P. integrifolia and P. axillaris, there were only 4539 homozygous SNPs, which constituted only 6.2 % of all polymorphisms, indicating that the self-incompatible P. integrifolia is highly heterozygous. Of all SNP loci, 60,701 were located within the CDS, of which 60,561 were biallelic; 9537 were located within 5ʹUTRs, of which 9512 were biallelic; 17,983 were located within 3ʹUTRs, of which 17,923 were biallelic. There is a small group of 786 SNPs whose location could not be determined.
As with the unfiltered SNPs, SNP frequency was calculated by dividing the total length of the reference transcriptome (Table 2) by the total number of SNPs (Table 5). Overall, the SNP frequency was 1/597 bp. The SNP frequency was 1/2056 bp between P. axillaris and P. exserta, 1/726 bp between P. axillaris and P. integrifolia, and 1/669 bp between P. integrifolia and P. exserta.
All SNPs among the three species were distributed across 20,606 unigenes (62.5 % of the total unigenes), corresponding to~76 % (40,556,099/53,135,953 bp) of the entire unigene length. When only considering the unigenes containing SNPs, the overall SNP frequency was 1 SNP/456 bp (89,007 SNPs/40,556,099 bp) and the highest SNP frequency was 1/89 bp. Among the SNPcontaining unigenes, only 1944 (9.4 %) had 10 or more SNPs (Table 5). This suggests that the SNPs were broadly distributed across the transcriptome, which might facilitate SNP marker selection for genome wide association (GWAS) studies.
SNPs between P. axillaris and P. exserta were distributed across 12,060 unigenes (Table 5). There were only 110 unigenes with 10 or more SNPs. The maximum SNP frequency per unigene was 1 SNP/111 bp. Between P. axillaris and P. integrifolia, SNPs were identified in 17,949 unigenes. There were 1466 unigenes with 10 or more SNPs. The maximum SNP frequency per unigene was 1 SNP/96 bp. There were 18,032 unigenes with SNPs between P. integrifolia and P. exserta, and 1722 unigenes with more than 10 SNPs. The maximum SNP frequency per unigene was 1/89 bp.
Due to the stringent nature of the SNP discovery parameters employed (i.e., filtering out three SNPs occurring within 100 bp of each other and filtering out SNPs within the first and last of 30 bp of a transcript, etc.), the above SNP frequencies calculated after the filtering steps are likely underestimated. Additionally, for organisms without a fully sequenced genome, using a transcriptome de novo assembly followed by variant detection can result in underestimation of expressed variants [52]. Even with the collection of different tissue types in our study, there are still likely to be undetected variants.
The percentage of unigenes with ten or more SNPs was 0.9 %, 8.2 % and 7.5 % of the entire SNP-containing unigene set between P. axillaris and P. exserta, P. axillaris and P. integrifolia, and P. integrifolia and P. exserta, respectively (Fig. 5). On average, there were 2.7 SNPs/ transcript, with 2.2, 0.78 and 2.4 SNPs/transcript between P. axillaris and P. integrifolia, P. axillaris and P. exserta, and P. integrifolia and P. exserta, respectively. When evaluating the transcripts with more than ten SNPs for GO annotation, the gene percentage in each category was in proportion with the overall transcriptome GO annotation (Additional file 1: Figure S5), indicating no particular set of genes were enriched regarding SNP frequency.
To gain additional insight into base substitutions occurring at each polymorphic site, the transition to transversion ratios (Ts/Tv) were determined (Table 6). Overall, the transition to transversion ratio (Ts/Tv) was 1.7:1. Across all species comparisons, the Ts/Tv was consistently highest for SNPs within CDS, followed by SNPs located in 3ʹUTRs, while the Ts/Tv was always lowest for SNPs located in 5ʹUTR. The overall Ts/Tv was relatively stable between species, with a range of 1.6 to 1.8 ( Table 6). The Ts/Tv ratio was used as a critical metric for assessing the specificity of new SNP calls in human genome research [53], and might be a useful parameter for subsequent Petunia SNP discovery. Normally, assuming that mutations are completely random, the Ts/Tv would be 0.5. Our Ts/Tv data indicate that each type of transitional change is produced more than three times as often as each type of transversion. Transitions occurring more frequently than transversions in transcriptome-derived SNPs have also been reported in other plant species, including Chinese pine [54] and melon [55]. The Ts/Tv bias could be the result of a high degree of methyl C to U in genomes [56,57]. In plants, Ts/Tv can vary across species. For example, in the exome assembly of four Neotropical tree species, Ts/Tv varied between 1.5 and 1.7 [58]. In the transcriptome from Norway spruce, the Ts/Tv was around 1.2 to 1.5 depending on the sequencing quality cut-offs [59]. It has been suggested that a higher C↔T transition can be accompanied by a higher number of Ts. The same situation was observed in our study. For instance, SNPs between P. integrifolia and P. axillaris, and between P. exserta and P. integrifolia, had a higher rate of C/T mutations than was observed between P. axillaris and P. exserta, which was accompanied by a slightly higher Ts/Tv. Minor allele read count frequency was used to calculate the frequency of short reads aligned to the least common allele in all three genotypes. Using a small number of genotypes for SNP discovery, this measure should provide better support for SNP confidence than using minor allele frequency (MAF), where the least prevalent allele frequency is calculated based on the genotype of each individual in a given population. For example, if there are only a few individuals in the panel, the smallest minor allele frequency will be 1/2n (n is the number of individuals), in our case, the minimum MAF would be 16.7 %. However, the average depth of coverage was above 100, and for RNA-seq data, the preferentially expressed genes will have a much higher read depth. Thus, even employing a depth cut-off value of 10, a heterozygous genotype call could have one allele with a much higher read count than the other allele, which might actually be caused by sequencing error or misalignment. This can be detected using the minor allele read count frequency. SNPs with minor allele reads count frequency ranging from 0.3 to 0.5 accounted for 61 % of the total SNPs (Fig. 6).

SNP validation
A total of 55 primer pairs were selected for PCR amplification. Of these, 50 amplified single DNA fragments, containing a total of 51 SNP loci, matching the expected sizes in P. axillaris and were selected for further sequencing (Additional file 1: Table S1). The failed/unclear amplification might be from primers located across intron-exon boundaries, or from the amplification being across a large intron or in paralogous genes. Three of the loci were monomorphic, five yielded poor sequence

CAPS marker design and validation
Twenty restriction enzymes were used to computationally predict CAPS markers from 15,701 SNP loci (Additional file 5: Dataset S4). Some SNP loci are predicted to be digested by two different restriction enzymes, with one recognizing the reference allele and the other recognizing the alternative allele, or by different restriction enzymes with overlapping recognition sequences. For example, for the restriction enzymes AluI and SacI, SacI is a 6 bp cutter recognizing/cutting sequences of GAGCTˇC and AluI is a 4 bp cutter that recognizes/cuts AGˇCT, thus all AluI recognition cites are also recognized by SacI. These analyses led to the identification of a total of 17,635 predicted CAPS markers (Supplemental dataset 4). Fourteen putative CAPS markers were randomly selected for validation.
All primer pairs amplified products of the expected size, and 11 were polymorphic. All 14 amplicons for P. axillaris and P. exserta had the expected genotype, while four amplicons from P. integrifolia were homozygous while the genotype prediction was heterozygous. However, the CAPS markers generated from this study have some limitations. First, the CAPS marker design pipeline only considered possible restriction enzyme recognition differences between species, the number/size of the fragments generated from digestion was not taken into consideration. Thus, there might be multiple small fragments yielded from the digestion, or the digested fragments between genotypes may only vary by a few base pairs. Either situation could not be scored successfully on agarose gel genotyping systems. Second, there could be sequence differences between the P. axillaris reference genome and the actual parental genotype used for this study (a different P. axillaris accession), and some restriction enzyme recognition sites predicted in the pipeline might not exist when screening our samples.
A total of 2949 primer pairs were developed from 3027 unigenes (Additional file 6: Dataset S5). The designed SSR primer pairs included 1,481 di-, 1,867 tri-, 88  Hexa-nucleotide tetra-, 16 penta-, and 47 hexa-nucleotide repeats ( Table 7). Forty-eight SSR markers were randomly selected for their utility in petunia (Additional file 7: Dataset S6). Of these, 36 (75 %) were successfully amplified and displayed polymorphisms between at least two of the three species examined. For example, 24 SSRs were polymorphic between P. integrifolia and P. axillaris, 28 were polymorphic between P. axillaris and P. exserta, and 24 were polymorphic between P. integrifolia and P. exserta. This success rate was higher than for our previous results developing SSR markers from P. axillaris ESTs [66].

Representation and in silico expression analysis of Petunia plastochron-related transcripts
Several genes from Arabidopsis are known to influence plastochron (Table 9). Local tBLASTN searches of the Petunia spp. transcriptomes using 14 Arabidopsis plastochron-related proteins as the query sequences revealed the presence of highly homologous unigenes within each transcriptome ( Table 9). The majority (34 out of 46) of the recovered transcripts are predicted to encode full-length proteins despite the fact that many plastochron-related genes, including ERECTA, ERECTA LIKE-1, SERRATE and SLOW MOTION, yield transcripts close to 3 kb or greater. Based on predicted amino sequence alignments, an additional five unigenes appear to be truncated by approximately 60 nucleotides (20 amino acids) or less at their N-termini (Table 9). This targeted analysis supports the transcriptome wide assessment of quality using CEGMA analysis (Fig. 1c) and indicates an overall high quality assembly of each transcriptome. Several of the previously characterized plastochron-related genes are preferentially expressed within the shoot apical meristem or developing leaf primordia [16,18,20]. Congruent with these previous findings, among the plastochron-related homologs identified in petunia, those related to AMP1, ERECTA, PIN1, TEL1, KLUH, SERRATE and SPL15 display enriched expression within shoot apices compared to the additional tissues examined (Table 9). These data suggest conservation of both gene content and expression pattern of plastochron-related genes between petunia and Arabidopsis.

Mapping plastochron-related homologous transcripts in an interspecific Petunia population
We identified Petunia transcripts homologous to numerous plastochron-related genes from Arabidopsis (Table 9). With the previously identified SNPs, 13 of these transcripts were converted to CAPS markers and used for linkage mapping (Additional file 1: Table S2). Together with the previously published SSR and CAPS markers [11], we constructed a linkage map with a total length of 289 cM consisting of 90 markers. Of those, eight were CAPS markers developed from plastochron-related transcripts located on five chromosomes (1, 2, 3, 6, and 7). The average linkage group length was 41.3 cM with a range from 32.7 cM (Chr7) to 58.7 cM (Chr5). The average marker density was 3.2 cM. Similar to our previous study [11], three QTL for development rate located on chromosomes 1, 2 and 5 were detected. Together they explained 37 % of the variation for development rate (Fig. 7; Table 10). However, the QTL location on chromosome 1 shifted away from a CAPS marker developed from the isopentenyl transferase gene SHO, a gene originally identified in an activation-tagged line exhibiting increased lateral shoot production [67]. We found that although markers developed from plastochron-related genes were located on Table 8 Summary of SSR repeat motif types and their corresponding repeat unit numbers for di-and tri-nucleotide repeats   Repeats  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  Total   AC/GT  -203  132  76  56  38  15  15  6  7  7  4  2  2  1  5 6 4   AG/CT  -209  130  81  40  34  25  17  17  14  13  10  4  5  9  3  6 1 1   AT/TA  -122  64  29  19  10  9  8  7  3  2   the same chromosomes as two of the QTL, they did not co-localize with development rate QTL regions. Linkage maps with higher resolution will help better understand the mechanisms for this trait.

Conclusions
High quality transcriptome assemblies were generated by RNA-seq data from three Petunia spp. and these were mined for molecular markers and plastochronrelated transcripts. We used a SNP discovery pipeline to identify over 80,000 SNPs across three species. As the majority of the SNPs were located in exons and over 80 % of the SNP-containing transcripts can be annotated, these SNPs could possess utility in future gene discovery and breeding applications. Thirteen percent of the identified SNPs were polymorphic among all three species, and can be further used for comparative genomics and population diversity studies. Additionally, over 10,000 SNPs were transformed into CAPS markers and an additional set of SSR markers were developed. Plastochron-related transcripts were identified, converted to CAPS markers, and genetically mapped in a P. integrifolia × P. axillaris F 2 population. Overall, these data provide sequence, expression information and a large marker resource for the Petunia community.

Plant material and growth conditions
Seeds of Petunia axillaris (PI667515) and Petunia exserta (OPGC943) were obtained from the Germplasm Resources Information Network. Seeds of Petunia integrifolia were purchased from Diane's Flower Seeds (Ogdon, UT). Five different tissue types were harvested for RNAseq: 1) whole seedlings, 2) callus, 3) shoot apices, 4) flowers of mixed developmental stages, and 5) trichomes. For whole seedlings, P. axillaris, P. integrifolia and P. exserta seeds were sown on 100 mm diameter Petri plates containing half-strength Murashige and Skoog (MS) media. Plates were left unsealed to allow air exchange and grown under a 16-h photoperiod provided by fluorescent lamps (ca. 75 μmol m −2 s −1 ) at 22°C. The date when seeds began to germinate was recorded (generally 2-3 d after sowing) and seedlings where harvested 21 d after this date. At harvest, media was removed from the roots, and plants were placed in 1.5 ml tubes and immediately frozen with liquid nitrogen. For callus generation, seeds of each species were sterilized by: 1) soaking in sterilized water for 30 min, 2) rinsing with 90 % ethanol, 3) soaking in 50 % bleach for 7 min and 4) rinsing with sterilized water. Twenty seeds per species were sown in a magenta box containing Coverage refers to whether a unigene encodes a predicted full-length protein (++), is truncated by~20 amino acids or less (+), or is truncated by greater than 20 amino acids (−) b The designation Paxi, Pexs and Pint refer to unigenes derived from P. axillaris, P. exserta and P. integrifolia, respectively  For the collection of shoot apices, flowers and trichomes, seed was sown in 288-cell plug trays in a soilless media and placed under in a greenhouse under intermittent mist at 24°C. After emergence of the first true leaf, seedlings were removed from mist and placed in a glass-covered greenhouse at 20°C under a 16 h photoperiod under ambient light, supplemented with 60 μmol m −2 s −1 provided by high-pressure sodium lamps. Shoot apices were collected when plants had unfolded 6-8 true leaves. Small leaves were removed with forceps to minimize contamination with leaf petiole tissue. The shoot apices from ca. 100 plants of each species were excised, placed in 15 ml conical tubes, and immediately frozen in liquid nitrogen. Flowers were harvested when a full range of reproductive stages of development were present, ranging from flower buds ca. 5 mm in length through the beginning fruit development. Harvested flowers were placed in paper envelopes, and frozen in liquid nitrogen. In order to harvest trichomes, stems and petioles from flowering Petunia plants were cut into ca. 5 cm pieces and immediately frozen in liquid nitrogen. The frozen tissue pieces were held using chilled forceps and the trichomes were removed from the stems by gentle scraping with a frozen paint brush and were collected in a mortar containing liquid nitrogen. Trichomes were immediately ground and transferred to 2 ml crew cap tubes for storage at −80°C.

RNA extraction
Total RNA was extracted from each tissue sample using the Trizol® reagent (LifeTechnologies) following the manufacturer's instructions. Forty micrograms of RNA was treated for DNA contamination using RNase-free DNase set (Qiagen). DNase-treated RNA was purified using the RNeasy® MinElute Cleanup kit (Qiagen). RNA yield and quality were determined using gel electrophoresis, spectroscopy and the Agilent 2100 BioAnalyzer RNA 6000 Pico chip with RIN values ≥ 8.0 achieved (Agilent Technologies).

RNA-seq library construction, sequencing, and analysis
A TruSeq RNA Sample Preparation kit (Illumina) was used to construct the cDNA libraries for sequencing. Illumina barcodes allowed pooling of 7-8 libraries per lane. Sequencing was completed on the Illumina HiSeq 2500 next-generation sequencer to 100 nt paired-end at the Michigan State University Research Technology Support Facility (RTSF; http://rtsf.msu.edu/; East Lansing, MI). Raw read quality was assessed with FastQC (http:// www.bioinformatics.babraham.ac.uk/ projects/fastqc/). Sequences were filtered and trimmed based on quality metrics and adapter sequences were removed with TrimmomaticPE [68]. The TrimmomaticPE options employed included ILLUMINACLIP, SLIDINGWIN-DOW:5:20, LEADING:5, TRAILING:5, and HEAD-CROP:10-14. TrimmomaticPE outputs sequences into paired-end and single-end files requiring no further mate-specific filtering. Cleaned reads were reassessed with FastQC for quality visualization to ensure no further filtering was required.

De novo assembly and expression analysis
Reads meeting quality standards were de novo assembled using the Velvet (version 1.1.07) and Oases (version 0.2.08) packages [69,70]. K-mer lengths of 55,57,59,61,63,65,67, and 69 were tested and metrics collected for each assembly. Criterion including the total number of transcripts produced, Velvet N50 length, average transcript length, and average transcript read coverage were considered. A k-mer of 61 was selected representing the best balance of metrics for all assemblies. Two de novo assemblies were completed for each Petunia transcriptome to reduce redundancy and limit the number of raw reads confounding assemblies. Reads from the callus, shoot apex and seedling libraries composed the first de novo assembly. Resulting transcripts were concatenated into a single pseudomolecule or 'artificial chromosome one' to serve as our core tissue reference transcriptome. Reads from the trichome and flower libraries were then  [46]. Reads aligning to the artificial chromosome were discarded to eliminate redundancy of data already assembled. All unmapped reads from the flower and trichome libraries were then combined with the original callus, shoot apex, and seedling reads in a second de novo transcriptome assembly [42]; Supplemental Fig. 1). Transcripts from the second and final assembly were filtered for low complexity. Representative transcripts, comprised of the longest isoform were extracted from the final assemblies. Completeness of the assemblies was assessed by mapping all reads from each tissue library back to their respective de novo assembly (per species) individually with Bowtie and TopHat [46]. CEGMA or Core Eukaryotic Genes Mapping Approach [43], a pipeline used to accurately annotate core genes in eukaryotic genomes, was also searched to determine the quality and completeness of the assemblies. exserta representative transcripts were removed as suspected contamination based on their annotation. Functional annotation of the transcripts was also queried using web-based FastAnnotator [73]. FastAnnotator facilitates the integration of the annotation tools Blast2GO, PRIAM, and RPS BLAST, providing gene ontology classifications (GO), enzyme and domain identification. Transcripts were translated in batch using OrfPredictor [74] and orthologous and paralogous proteins in the three species were assigned with OrthoMCL v2.0.2 with default parameters and an e-value cut-off of 1e-10 [50].

SNP Identification
Sequence reads from each tissue type were pooled for each species and mapped to the Petunia axillaris transcriptome using the Burrows-Wheeler Aligner (BWA) program [75] with default values. Duplicate reads were removed after the initial alignment, to eliminate reads that mapped to the same position of the transcriptome. Duplicate removal was performed for the aligned reads using picardTools/1.89 (broadinstitute.github.io/picard). The subsequent local realignment to correct misalignments due to the presence of indels was performed using the Genome Analysis Toolkit (GATK) [51] with the IndelRealigner function. Bam files for each genotype were compressed using ReduceReads from GATK. The initial variant calling was performed by HaplotypeCaller from GATK with a Phred-scaled confidence threshold of 30. After the initial variant calling, SNPs were filtered based on several criteria. The first round of SNP filtering employed the following steps: exclude three SNPs within 100 bp; filter out variants with zero mapping quality constituting more than 10 % of all reads at that site; exclude SNPs covered only by sequences on the same strand with an FS value (Phred-scaled p-value using Fisher's exact test to detect strand bias) >60; exclude SNPs with a minor allele frequency <0.01; and exclude SNPs with low QD (quality by depth), low quality (<11), and low total depth of coverage (<11) by the default parameters recommended from GATK best practice with the VariantFiltration tool [53,76]. The second round of SNP filtering was performed by custom Python (Python 2.7.8) [77] scripts with the following criteria: SNPs within the beginning and end 30 bp of the reference transcripts were excluded; SNPs were selected with at least two genotypes having a depth of coverage greater than 10; exclude loci with a heterozygous genotype call in Petunia axillaris. The intron-exon boundaries were predicted by aligning the Petunia axillaris transcripts to the draft Petunia axillaris genome scaffold sequences with GMAP [78]. SNPs within 30 bp of the exon-intron boundary region were then filtered out. Indels were not called because alternative splicing may impede reliable indel discovery. A custom Python script was used to extract the 5ʹ-UTR and 3ʹ-UTR and to calculate the distribution of SNPs in these regions. Depth of coverage was calculated by BEDTools [79]. Annotated unigenes containing SNPs were visualized by BGI WEGO (web gene ontology annotation plotting) (http://wego.genomics.org.cn/cgi-bin/wego/index.pl). KEGG pathways were assigned to unigenes containing SNPs using the online KEGG Automatic Annotation Server (KAAS) (http://www.genome.jp/tools/kaas/) [80]. KEGG Orthology (KO) assignment was applied using Bi-directional Best Hit (BBH) method with plant organisms (Arabidopsis thaliana (thale cress), Arabidopsis lyrata (lyrate rockcress), Theobroma cacao (cacao), Glycine max (soybean), Fragaria vesca (woodland strawberry), Vitis vinifera (wine grape), Solanum lycopersicum (tomato), and Oryza sativa L. ssp. japonica (Japanese rice) (RefSeq)) as references.

SNP validation
A total of 55 SNPs were selected for validation. Primers were designed by first choosing SNPs where at least two species exhibited polymorphism. Unigenes with exons greater than 600 bp were then selected, SNPs or alleles in the sequences were masked using IUBI./IUPAC nucleotide acid code. Batch Primer3 (http://probes.pw.us da.gov/cgi-bin/batchprimer3/batchprimer3.cgi) was used for primer design using the "SNP flanking primers design" option. The expected amplicon sizes ranged from 250 to 350 bp, with primer sizes ranging from 18 to 27 bp, and primer T m ranging from 57 to 63°C. Genomic DNA was extracted using the DNeasy plant mini kit (QIAGEN, Valencia, CA, USA). PCR amplification of genomic DNA was carried out in a 10 μl reaction containing 1 × PCR buffer, 0.2 μl 10 μM dNTPs, 0.6 U of DNA polymerase, and 5 ng template DNA.

CAPS and SSR marker design
SNP markers were transformed into CAPS markers using the following pipeline. First, the SNP output was transformed into VCF (variant call format) format. Then, by using the previous GMAP output as reference, a new VCF file was generated by transforming the current SNP locations and genotypes to corresponding SNP locations and variant calls in the P. axillaris draft genome scaffolds by a custom Python script. Primers were designed with the galaxy-pcr-markers pipeline (https://github. com/cfljam/galaxy-pcr-markers) with the modification of producing only one pair of primers, minimum amplicon size of 200 bp and maximum amplicon size of 400 bp. A total of 20 commonly used restriction enzymes: AluI, ApaI, BamHI, BbrPI, BfrI, ClaI, DdeI, DpnII, DraI, EcoRI, HaeIII, HincII, HinfI, HpaI, PvuII, RsaI, SacI, Sau3AI, SmaI, and TaqI were selected for the pipeline. The unigene dataset were used for SSR identification with MISA (microsatellite identification tool) (http://pgrc.ipk-gatersleben.de/misa/). The SSR identification criteria were 6, 5, 5, 5, 5 repeats for di-, tri-, tetra-, penta-, and hexanucleotides, respectively. Primer pairs were designed from primer3_ core (Primer3 2.3.6) (http://primer3.sourceforge.net/releases.php). Primer parameters were: designated amplicon size 100-280 bp, primer annealing temperature 55 to 65°C, primer length 18 to 28 bp, and GC content from 45 to 55 %. Unigenes with homology to known genes involved in specifying plant development rate [16-18, 20, 25] were searched for possible CAPS markers from the entire CAPS output. The CAPS markers were used to genotype the three species by following the previous PCR protocol. CAPS markers were digested by the above mentioned restriction enzymes (New England Biolabs, Beverly, MA) and separated on 2 % NuSieve GTG agarose (Lonza, Basel, Switzerland) plus 1 % UltraPure agarose (Life Technologies, Carlsbad, CA) gel with 100 V for 2 h at room temperature.
Candidate gene mapping and QTL mapping of development rate in an interspecific hybrid Petunia population An interspecific hybrid F 2 population developed from a cross between P. integrifolia and P. axillaris containing 164 individuals was used to identify QTL for development rate. Population development, the measurement of development rate, and the SSR-based genetic linkage map were reported previously [11]. CAPS markers developed from P. axillaris plastochron-related transcripts based on SNPs between P. integrifolia and P. axillaris were used to screen the F 2 population. If no CAPS were previously designed from our SNP set, CAPS markers were manually designed with CAPS Designer (http:// solgenomics.net/tools/caps_designer/caps_input.pl) and Primer3Plus (http://www.bioinformatics.nl/cgi-bin/pri mer3plus/primer3plus.cgi/). The updated genetic linkage map was generated with JoinMap 4.0 [81] with the Kosambi mapping function [82]. The recombination threshold was set to 0.3 and the logarithmic odds (LOD) score was set to 5. The locus order was calculated by the regression module of JoinMap4. Linkage group numbers were assigned according to the previous study [11]. QTL mapping was conducted with MapQTL 6 [83]. A permutation test with 1000 permutations was used to establish the LOD significance threshold. The initial QTL mapping was performed by interval mapping (a single-QTL model). Multiple QTL model mapping (MQM) was then used to reduce the background