Chloroplast Genome

The reference-guided assembly of the chloroplast genome was completed in alignreads using the chloroplast genome sequence of oleander (Nerium oleander: Apocynaceae; M. Moore and D. Soltis, unpublished data; see [17] for coding sequences) with only one copy of the inverted repeat included as the reference and no masking. A second alignment of the YASRA contigs was completed using Mulan [35], and the resulting consensus sequence was compared to the alignreads output and tested as an alternate reference input for alignreads. All of the consensus sequences were aligned to each other and the oleander reference using MAFFT 6.240 [36]. A preliminary annotation of the A. syriaca chloroplast genome was completed using DOGMA [37]. The annotations of protein coding regions, ribosomal DNA (rDNA), and sequences corresponding to transfer RNAs (tRNAs) were improved by fixing annotation and assembly mistakes (e.g., no start or stop codon, internal stop codons, frameshifts, endpoints of tRNA sequences) through examination of the actual assembly of the short reads at that location in Tablet vers. 1.10.09.20 [38] and/or comparisons to ortholog sequences from other asterids obtained from NCBI GenBank (e.g., Nerium, Coffea, Solanum). Various versions of the chloroplast genome sequence were used as a reference in alignreads to further refine the sequence based on new assemblies. Discrepancies between the assembly consensus sequences in introns and intergenic spacers were explored at the microread-level in Tablet to identify the correct sequence. An Asclepias tuberosa sequence [GenBank:GQ248251.1] served as a reference in alignreads to correctly assemble the trnH - psbA region.

The final annotation of the A. syriaca chloroplast genome was accomplished using BioEdit 7.0.5.3 [40] and visualized using GenomeVx [41]. The repeat content of the chloroplast genome was determined using the NCBI BLAST tool blast2seq. Only repeat motifs of greater than 30 bp and 85% sequence identity were characterized. Sequence identities were calculated in BioEdit on sequence alignments produced in MAFFT with the gapped positions removed. Ka/Ks values for pseudogenes were calculated using the Bergen Center for Computational Science Ka/Ks calculation tool http://services.cbu.uib.no/tools/kaks.

Mitochondrial Genome

The mitochondrial genome of A. syriaca was assembled using the alignreads pipeline and the 430,597 bp mitochondrial genome sequence of Nicotiana tabacum [GenBank:BA000042.1] [42] as a reference. Mitofy [43] was used to characterize the complement of protein coding, rRNA, and tRNA genes in the resulting contigs. The annotation of the mitochondrial genome of N. tabacum was updated using Mitofy before comparison of the gene content of the two genomes. BLAT [44] under default parameters was used to estimate the percentage of each A. syriaca mitochondrial gene assembled by calculating the proportion of the total length of the N. tabacum mitochondrial genes with hits in the A. syriaca contigs. The presence or absence of the cox1 group I intron in Asclepias was determined by using the Hoya sikkimensis (Apocynaceae) cox1 sequence [GenBank:AJ247588.1] as a BLAT reference.

Ribosomal DNA

The nuclear rDNA cistron, including 18S, 5.8S, and 26S, was assembled in alignreads using a 669 bp partial sequence from A. syriaca [GenBank:AM396892] spanning the internal transcribed spacers (ITS) as the reference. Following the initial assembly, the assembly was repeated using the reverse complement of the reference sequence as the new reference to correct for an inconsistency in YASRA that caused the assembly result to be strand specific. A consensus sequence of the contigs obtained through the first two analyses was then used as the reference for a final assembly. A 5S rDNA sequence for A. syriaca was assembled using a 330 bp sequence including the 5S locus from Solanum iopetalum [GenBank:AJ226045] as a reference in both the forward and reverse directions.

To assess intragenomic polymorphism, the complete genomic read pool was filtered for reads matching the consensus sequences of the assembled rDNA loci using BLAT and custom scripts [45]. Reads were discarded if their average Phred score was less than 20, and the remaining reads filtered so that any base in a read with a Phred score less than 20 was converted to 'N' using the FASTX-Toolkit [46]. Modified reads were assembled against the reference sequence using BWA [47], and at each position the number of reads differing from the consensus was tallied and a proportion calculated (script available from KW). In order to estimate the error introduced by the sequencing method, the same procedure was applied to the PhiX174 control lane of the sequencing run.

Positions with 2% or more of reads differing from the consensus were considered to represent true polymorphism. This cutoff is comparable to error rates previously found using a similar quality-filtering scheme [48], but may be somewhat conservative as it is much higher than the error seen in the PhiX lane. See additional file 1: Polymorphism detected among Illumina reads for the PhiX sequencing control. All positions of the PhiX genome exhibited <0.7% differing reads except position 3132, which had >17% of reads differing from the consensus. While the reason for this anomaly is uncertain, the overall average error rate in the PhiX lane is <0.05%.

De novo assembly

The Chloroplast Genome

The chloroplast genome of A. syriaca is 158,598 bp [GenBank:JF433943], excluding two small, unresolved regions that were not able to be assembled or Sanger sequenced due to polynucleotide stretches (rps8-rpl14 intergenic spacer and ψycf1), and has an inverted repeat (IR) of 25,401 bp (Figure 1). The initial assembly produced using the alignreads pipeline and the oleander reference contained 51 contigs with a median read depth of 246× and N50 of 4,683 bp. The longest contig was 14,884 bp. The final assembly using the finished A. syriaca sequence as a reference contained 22 contigs, had an N50 of 9,030 bp, and longest contig of 28,186 bp. The use of chloroplast contigs from 80 bp reads from other A. syriaca individuals (S. Straub and A. Liston, unpublished data) in combination with Sanger sequencing of select regions in the 0.5× genome individual and these other individuals, resulted in the addition of 6,622 bp and deletion of 1,402 bp of sequence. That large insertions and deletions were found relative to the original assembly was expected due to the limited power of reference guided assembly algorithms to reconstruct these differences. Only 38 bp of sequence from the original reference guided assembly were determined to be incorrect, producing substitution errors in the genome sequence. The sum total of these changes resulted in the alteration of approximately 5% of the total sequence length, indicating that 95% of the chloroplast genome of A. syriaca was assembled at 0.5× average genome coverage with reads of only 40 bp and an oleander reference with 85% sequence identity (considering only one copy of the inverted repeat). Although the original assembly was very good overall, comparison of the final assembly with the finished reference sequence highlighted the limitations of using 40 bp vs. 80 bp reads for sequence assembly. Even with a correct reference sequence, some regions were still not assembled properly (e.g., accD, ycf1) and some assembly mistakes from the original assembly were recreated (e.g., an erroneous 5 bp deletion in the middle of rbcL causing a frameshift).

The A. syriaca ψclpP has several features that indicate that it is likely a pseudogene. The nucleotide and protein sequences are divergent from those of oleander (81.3% exon nucleotide and 67.0% protein sequence identity) and coffee (80.8% exon nucleotide and 67.5% protein sequence identity). In comparison, oleander and coffee have 97.2% exon nucleotide and 98.9% protein sequence identity. There is also a 28 bp insertion near the end of exon 3 of A. syriaca (validated by Sanger sequencing) that causes a frame shift and the introduction of a stop codon. Sequence comparisons showed that clpP in oleander and coffee are likely under negative selection (Ka/Ks = 0.07 and 0.03 respectively); however, the Ka/Ks for Asclepias is 1.15 indicating that this locus is likely evolving neutrally. While the clpP introns have been lost multiple times in angiosperms [e.g., [63–65, 67, 68]], the coding region has been lost much less often. Other examples of loss of the coding region are only known from Scaevola (Goodeniaceae), Passiflora (Passifloraceae), Trachelium (Campanulaceae), Geranium and Monsonia (Geraniaceae) [67–69].

The accD gene of A. syriaca is also likely a pseudogene. This gene is 2,584 bp in A. syriaca compared to 1,497 bp in oleander. There is a repeat-rich (Table 5), approximately 1 kb insertion in the middle of the gene. Although the insertion does not interrupt the reading frame, the remainder of the sequence that is alignable with oleander is highly divergent. This gene can be highly variable, even within a genus [66], but due to the extent of sequence divergence, and lengths of the insertion and repeat motifs, this gene is likely non-functional in Asclepias. However, in pepper, a smaller, repeat-containing insertion of 144 bp in accD does not cause a frameshift nor prevent transcription of the gene [70]. There are numerous other examples of loss of accD from the plastome among angiosperms [e.g., [64, 65, 67–69, 71]].

The assembly of ψaccD proved especially challenging for Asclepias due to sequence divergence and the large insertion relative to the oleander reference. De novo assembly of 80 bp reads from another A. syriaca individual and Sanger sequencing confirmed the correct sequence. Further inspection indicated that the assembly errors in this region arose due to the presence of an accD pseudogene in the mitochondrial genome of A. syriaca, also confirmed by Sanger sequencing, which was more similar in sequence and length to the oleander chloroplast accD. A similar problem involving accD and a mitochondrial genome sequence was encountered in the assembly of the pepper plastome [70]. Problems of this sort are likely to arise for assembly of chloroplast genomes, especially with short reads, for regions of the genome that have been duplicated in the mitochondrial or nuclear genomes [72].

A third possible pseudogene in the A. syriaca chloroplast genome is ψycf1. This gene contains one of the two polynucleotide repeats of unresolved length in the genome and has a highly divergent sequence compared to oleander and coffee. The sequence of this gene was also checked and improved using de novo contigs and Sanger sequencing to eliminate the possibility of assembly errors as a cause of sequence divergence. Although ycf1 is one of the most quickly evolving chloroplast genes in seed plants [e.g., [66, 73, 74]], the full length copy of ycf1 in Asclepias is so divergent that it may be non-functional. The sequence identity between ycf1 in oleander and coffee is 87.5%, whereas the sequence identities between each of those species and Asclepias are only 79.4% and 72.2% respectively. In addition to decreased sequence identity, there were 62 and 61 indels inferred from alignments of the Asclepias ycf1 sequence with the oleander and coffee sequences respectively, whereas there were only 16 indels inferred from an alignment of oleander and coffee. Among angiosperms ycf1 has also been lost multiple times [e.g., [65, 67–69, 71]].

All three of the putative pseudogenes detected in the A. syriaca plastome, ψclpP, ψaccD, and ψycf1 have functions that are known to be essential for plant development and survival [75–77]. Additional studies to evaluate the functionality of these putative pseudogenes will be needed to confirm their designation as such, in combination with further analysis of the nuclear genome to locate functional copies of these indispensible genes or determine which other nuclear genes might encode a functional replacement as in, for example, the case of accD in grasses [e.g., [78, 79]].

The Mitochondrial Genome

A total of 130,764 bp of the A. syriaca mitochondrial genome was assembled in 115 contigs with a median read depth of 32× and N50 of 1,666 bp. See additional file 4: Asclepias syriaca mitochondrial genome contigs (also accessible at http://milkweedgenome.org). The longest contig was 6,948 bp. The extent to which the portion of the genome assembled for A. syriaca represents the whole mitochondrial genome is unknown due to the extensive variation observed for genome size, gene content, and organization among angiosperms [43, 80, 81]. Divergence from the organization of the Nicotiana reference likely prevented a more complete and contiguous assembly. Complete (rrn5) or nearly complete (rrn26, rrn18) ribosomal DNA sequences were identified among the contigs. Twenty-eight putative tRNA sequences representing 25 different tRNAs were also identified. Absence of a complete complement of tRNAs was not unexpected because the number of tRNAs in plant mitochondrial genomes is variable due to the replacement of native mitochondrial tRNAs with versions of chloroplast origin or import of nuclear encoded tRNAs of chloroplast or nuclear origin [80, 81]. For example, among the tRNAs in the A. syriaca mitochondrial genome, both a mitochondrial version and a chloroplast version of trnE- TTC were detected.

The protein coding gene content observed for A. syriaca was compared to that of the mitochondrial genome of Nicotiana tabacum because it is the only other asterid with a sequenced mitochondrial genome. The assembled part of the genome contained the same 37 protein coding genes of known function as the N. tabacum mitochondrial genome with an average coverage of 94.8%. Putative pseudogenes of rps2, rps7, and rps14 were identified in the A. syriaca mitochondrial genome contigs. Both the mitochondrial genome of N. tabacum and A. syriaca appear to lack functional copies of these three genes, in addition to lacking rpl6, rps8, and rps11, the genes lost in all angiosperms, but present in Marchantia[42]. Evidence of ψrps14 is present in tobacco, but there is no remaining evidence of the other two genes. Due to the pseudogenization or absence of these genes in the N. tabacum reference used for the A. syriaca assembly, a more complete characterization of the mitochondrial genome based on additional data will be required to determine if these genes are truly non-functional in A. syriaca or if their absence is an assembly artifact.

Extensive horizontal gene transfer has been hypothesized for the group I intron present in the cox1 gene of angiosperm mitochondrial genomes, with six horizontal transfer events hypothesized in Apocynaceae alone [82]. Alternatively, it has been argued that the distribution of the cox1 intron in plants more likely involved a single horizontal transfer followed by intron loss in many lineages [83]. All but one of the Apocynaceae genera sampled have contained the cox1 intron [82, 83], and it was also found to be present in the assembled cox1 gene of A. syriaca. This result adds to the evidence supporting the intron loss hypothesis over horizontal gene transfer in Apocynaceae.

Repetitive Element Characterization

A total of 564 repetitive elements were detected in the A. syriaca de novo contigs (Figure 2). The majority of these repeats were Ty1/copia- like and Ty3/gypsy-like LTR retrotransposons (Figure 2), which are abundant in plants and common across eukaryotes [84–86]. Only the most highly represented repeats in the genome would have high enough coverage at 0.4× total nuclear genome coverage to be assembled, so it is not surprising that the majority of the repeats were Ty1/copia- like and Ty3/gypsy-like retrotransposons. In plant species where retrotransposons have been characterized, it appears to be lineage specific whether Ty1/copia-like or Ty3/gypsy-like retroelements are more abundant in the genome, likely because retrotransposons can proliferate or be lost from genomes over relatively short evolutionary time scales [85, 87]. Among other asterids, Ty3/gypsy-like retrotransposons are more prevalent in the sunflower (Helianthus annuus), tomato, and potato genomes, while in carrot (Daucus carota), as in A. syriaca, Ty1/copia-like retrotransposons are more numerous [88–90].

The numbers of retroelements reported here are certain to be different than the actual numbers of retroelements in the sampled A. syriaca genome because the de novo contigs were relatively short compared to the length of retrotransposons, possibly resulting in different regions of the same retroelement being counted multiple times. Additionally, six contigs had more than one good repeat type match, slightly inflating the numbers discovered. Only approximately 2% of the putatively nuclear reads had hits to the repeat-containing contigs. The estimate of raw reads with hits to the plant repeat database did not give a better estimate of the total repeat content of the A. syriaca genome than the hits to assembled repeats with only 0.24% of reads with one or more hits. These results indicate that at 0.4× genome coverage, only an initial characterization of the most common repeats present in a plant genome can be made, while an estimation of the overall repeat content of the genome requires deeper coverage.

To be conservative in selecting 25 candidates for testing, loci from the 56 contigs with putatively mitochondrial or chloroplast sequence were removed from the pool of candidates, leaving 209 primer sets for putatively nuclear microsatellite loci. See additional file 2: Primer sets designed using BatchPrimer3 for 184 nuclear microsatellite loci in Asclepias syriaca not tested for amplification success. Of the 25 primer sets tested, all successfully amplified in A. syriaca (Table 6). It is important to note that with the low overall genome coverage, at least some of these microsatellite-containing contigs are likely to have originated from multiple, nearly-identical paralogous loci. As is the case with all microsatellites, the ultimate proof in the utility of these markers for identifying intraspecific variation (e.g., linkage mapping or population genetics) will come from obeying a Mendelian locus/allele model that meets expected transmission ratios (in the case of known pedigrees) or Hardy-Weinberg expectations (in the case of wild populations). If the markers show evidence of fixed heterozygosity or dominant-only variation, this would constitute evidence that our microsatellite 'contigs' are the product of microreads representing two or more paralogs. These loci are currently being examined in greater detail to determine their locus specificity.

The source of the rDNA repeats identified by RepeatMasker in the de novo contigs (Figure 2) that remained after the removal of A. syriaca rDNA microreads was also explored. BLAST searches (blastn) of the NCBI nucleotide collection revealed that their source was likely plant fungal pathogens or endophytes. For example, the top BLAST hits for two of theses contigs were to an Alternaria alternata gene for large subunit (LSU) rRNA (1e^-105) [GenBank:AB566324.1] and a Davidiella tassiana partial LSU rRNA gene sequence (1e^-92) [GenBank:FN868880.1]. Approximately 0.02% of the total reads could be categorized as fungal rDNA when the sequences for the top GenBank hits were used as a BLAT database. Detection of fungal contaminants of a plant genomic DNA preparation is unsurprising due to the ubiquitous association of these two types of organisms.

Ribosomal DNA

A 385 bp contig containing a 120 bp 5S rDNA sequence [GenBank:JF312047] and an rDNA cistron of 7,541 bp [GenBank:JF312046] were assembled for A. syriaca (Figure 3). Repeat structure in the external transcribed spacer (ETS) and intergenic spacer (IGS) regions of the rDNA cistron prevented further extension of that contig and produced an error in the assembly where reads corresponding to similar repeats were incorrectly piled up, highlighting one of the primary difficulties of genome assembly using short reads [32, 91]. Consequently, the first 280 bp of the contig sequence were removed prior to downstream analyses. The median read depth for the final assemblies were 406× for 5S and 738× for the 18S-5.8S-26S cistron. Using the median read depth as an estimate of the number of rDNA copies sequenced and the nuclear genome coverage of approximately 0.4×, rough approximations of the number of 5S rDNA repeats and of rDNA cistron copies are 1,015 and 1,845 respectively. These estimates are in line with values observed for other plants [92, 93].

Fully characterized rDNA cistrons in plants have at least one TATA box upstream of 18S, which demarcates the boundary between the ETS region and the repeat-rich non-transcribed spacer [94, 95]. However, repeat units may also appear between the transcription initiation site (TIS) and 18S [96]. The conserved plant TIS sequence was not observed in the A. syriaca sequence upstream of 18S, so the entire ETS region was not assembled due to internal repeated sequence. However, another ETS motif conserved in plants and always located approximately 420 - 450 bp upstream of the start of 18S [TGAGTGGTGG; [97]] was located 423 bp upstream of 18S in the sequenced portion of ETS. Due to the length variation in the ETS regions in other asterids (Solanum ~1000 bp; Nicotiana >3000 bp in N. tabacum; Olea ~2000 bp), it is difficult to estimate the completeness of our assembly for this region and the cistron as a whole, which is between 9 and 12 kb in Nicotiana[96, 98, 99].

For intragenomic rDNA polymorphism characterization, quality filtering by average quality (q = 20) removed zero reads from the 5S pool, and 148 reads from the 18S-5.8S-26S cistron pool (0.1%). Masking the remaining nucleotides with qualities below 20 affected 1.7% and 6.2% of nucleotides in the two pools, respectively. A high percentage of positions in the 5S locus were polymorphic (Figure 4). Most of these positions had less than 8% of reads differing from the consensus, though some ranged as high as 14%. Conversely, only 19 positions in the 18S-5.8S-26S cistron were polymorphic (Figure 4). Polymorphisms were found in both coding and spacer regions, although none were found in either the 5.8S or ITS2 regions. While spacer regions contained a slightly higher proportion of polymorphic positions than coding regions, positions with the greatest proportion of differing reads were found within coding regions, in particular in the 26S region.

Levels of intra-genomic polymorphism among copies of the rDNA loci generally agree with previously measured and theoretically expected levels. The high proportion of polymorphic positions in 5S is concordant with previous observations for this locus; however, the levels measured here are somewhat higher than those reported in other plant and animal species [e.g., [100–102]]. Intra-genomic polymorphism among copies of the rDNA cistron in this individual is comparable to rates found in previous studies characterizing polymorphism within Drosophila[103] and fungi [104]. While the 19 sites found in this study are within the previously observed range of variation, comparisons are difficult due to the difference in copy number between these taxa (45-180 in the fungi, 200-250 in Drosophila, approximately 1800 in this study), differences in the definition of polymorphic loci, and differences in sequencing technologies (whole-genome shotgun sequencing vs. Illumina).

The distribution of polymorphic positions matches theoretical expectations, with a higher proportion of positions found within spacer rather than coding regions. The very high proportion of differing reads found at two positions in the 26S region may be caused by their location within expansion segments of this region, expansion segments 9 and 8, respectively [105]. The lower levels of polymorphism observed for the 18S-5.8S-26S cistron relative to the 5S rDNA loci are concordant with the molecular evolutionary processes associated with each rDNA type. While strong concerted evolution and strong purifying selection act to homogenize the 18S-5.8S-26S cistron copies within and among repeat arrays, weak concerted evolution and birth-and-death processes for 5S rDNA allow accumulation of polymorphism [100, 106].

Protein Coding Nuclear Genes

From 19,899 Catharanthus ESTs, 9,266 unigenes were created. Of these, 72 corresponded to chloroplast, mitochondrial or rDNA genes and were removed from the data set. Among the remaining repeat masked unigenes, 8,077 had between one and 138 unique microread hits. The median number of hits per kilobase of sequence was 7.17 (95% CI 7.03 - 7.29). The median coverage of the unigenes was 0.29× (95% CI 0.28× - 0.29×; Figure 5a). The 618 genes with coverage levels above 0.984× were outliers. Many of these genes likely correspond to members of multigene families or other duplicated genes, increasing the chances of non-orthologous BLAT hits in conserved domains. For example, a unigene with 6.79× coverage had multiple BLAST hits (E = 0) to eudicot H⁺-ATPases (P_3A-type). Within the P-type ATPase gene super family, the P_3A-type subfamily has 10 and 11 members in rice and Arabidopsis respectively [107]. Based on the overall genome coverage and the unigene coverage, the A. syriaca genome likely contains a similar number of P_3A-type ATPases.

Another explanation for the high coverage of outlier loci was indicated by the masked unigene sequences themselves. Many of the unigenes with the highest coverage still contained sequences that could produce multiple non-homologous hits (e.g., internal polynucleotide stretches), which in turn produced stacks of reads. Read overlap due to this phenomenon inflated the coverage values, so many of the very high coverage values are likely artifacts of the unigene creation and masking process. The presence of PCR duplicate reads in the data set could also have caused read overlap and overestimation of coverage. Most of these reads should have been excluded by counting only unique BLAT hits, but any sequencing errors in the duplicates or length differences due to differing sequence quality, especially at the ends of the reads, could produce nearly identical hits.

Of the 1,086 COSII alignments 721 had between one and 23 unique A. syriaca microread hits with a median gene coverage of 0.14× (95% CI 0.13× - 0.15×; Figure 5b) and median hits per kilobase of sequence of 3.40 (95% CI 3.13 - 3.64). The 35 genes with greater than 0.473× coverage were outliers. There is a possibility that some of these genes are single copy in other asterids, but duplicated in Asclepias, leading to the higher coverage levels. More likely, as in the unigene analysis, sequence features, including microsatellites that were not masked in the COSII alignments, and PCR duplicates led to increased estimates of coverage. Of the five genes with greater than 1× coverage, the only BLAT hits for two were to microsatellites and a third had multiple microsatellite and a single other hit.

The median coverage values for the unigene and COSII data sets would be expected to be close to the estimate of 0.4× coverage for the nuclear genome. However, the observed values deviated from the expected coverage. One explanation for this is that the overall coverage of the A. syriaca genome was slightly overestimated due to the presence of PCR duplicates in the read pool. It is also possible that the genome size of the sequenced individual was larger than the genome sizes estimated from the experimental population, again causing an overestimate of genome coverage. Divergence from the reference sequences would also cause an underestimate of coverage for A. syriaca orthologs. The lower amount of divergence from Catharanthus, a member of the same family, is reflected in the higher estimated coverage for the unigenes than the COSII alignments, which contained sequences of more distantly related asterids. The coverage value for the unigenes may also be higher due to the reasons discussed above. Another feature of both reference data sets that would cause an underestimate of coverage was the absence of introns. The microread sequences from genomic DNA would also represent these features of the genes, but would not produce BLAT hits for these genes. Reads spanning intron-exon boundaries would also have been excluded in many cases due to the mismatch of the intron sequence with the following exon in the unigenes and COSII alignments. Even with these issues confounding coverage estimates, these two explorations of the nuclear gene content of the A. syriaca genome demonstrate that even at 0.4× genome coverage, much of the protein coding portion of the nuclear genome can be detected.

Even with low overall coverage of the 721 COSII genes with A. syriaca microread hits, the alignments were still useful for marker development for Asclepias and Apocynaceae molecular phylogenetics. The 28 genes that had 0.5× or greater coverage were evaluated for their suitability for this purpose. Microreads flanked putative intron locations in 15 genes of this set, allowing primer design for amplification of one or more introns per gene (Table 3).