Comparative assessment shows the reliability of chloroplast genome assembly using RNA-seq

Chloroplast genomes (cp genomes) are widely used in comparative genomics, population genetics, and phylogenetic studies. Obtaining chloroplast genomes from RNA-Seq data seems feasible due to the almost full transcription of cpDNA. However, the reliability of chloroplast genomes assembled from RNA-Seq instead of genomic DNA libraries remains to be thoroughly verified. In this study, we assembled chloroplast genomes for three Erysimum (Brassicaceae) species from three RNA-Seq replicas and from one genomic library of each species, using a streamlined bioinformatics protocol. We compared these assembled genomes, confirming that assembled cp genomes from RNA-Seq data were highly similar to each other and to those from genomic libraries in terms of overall structure, size, and composition. Although post-transcriptional modifications, such as RNA-editing, may introduce variations in the RNA-seq data, the assembly of cp genomes from RNA-seq appeared to be reliable. Moreover, RNA-Seq assembly was less sensitive to sources of error such as the recovery of nuclear plastid DNAs (NUPTs). Although some precautions should be taken when producing reference genomes in non-model plants, we conclude that assembling cp genomes from RNA-Seq data is a fast, accurate, and reliable strategy.

The development of high-throughput sequencing technologies has led to a rapid increase in the availability of chloroplast genomes 16,[22][23][24] making possible the use of complete molecules in phylogenomic analyses 17,[25][26][27] . At present, more than 2,500 complete chloroplast genomes are available 28 . However, the use of complete genome sequencing to obtain reliable chloroplast genomes also poses some caveats and remains relatively expensive. Transcriptome sequencing (RNA-Seq) is comparatively less complex because it yields only the sections of the genome that are transcribed into RNA, providing a relatively cheap and fast method to obtain large amounts of functional genomic data [29][30][31][32] . Accordingly, global initiatives such as the 1,000 plants (1KP) project have generated a wealth of transcriptomic data for over more than 1,000 plant species 33 . Since the chloroplast genome appears to be fully transcribed, RNA-Seq data could potentially be used to obtain the complete chloroplast genome 34 . However, the reliability of assembling chloroplast genomes from transcriptomic versus genomic data has not been thoroughly evaluated.
In this study, we compared the reliability of RNA-Seq to genomic DNA libraries to obtain cpDNA complete sequence. For this purpose, we assembled for the first time the complete chloroplast genome of three Erysimum (Brassicaceae) species: Erysimum mediohispanicum, E. nevadense, and E. baeticum from genomic libraries.
Erysimum constitutes an interesting case study because it is a genus that encompasses wide diversity attained through rapid and complex evolutionary processes [35][36][37] , while being evolutionarily close enough to Arabidopsis thaliana to render the use of genomic and transcriptomic references from this model species relatively easy. We assembled the chloroplast genomes of these three species using different computational approaches and compared several genetic features (gene content, presence of repeats, microsatellites -SSRs-, etc) across genomes obtained RNA-Seq or genomic DNA. Based on these results, we assessed a) The characteristics of the chloroplast genomes of Erysimum spp.; b) the genomic coverage provided by RNA-Seq across species and c) a bioinformatic approach to ensure reliable chloroplast genome assembly from transcriptomic data. In the light of these results, we propose a pipeline-like methodology for processing RNA-Seq reads into high-quality cp genomes.

Material and Methods
Plant materials. Fresh leaves and flower buds of Erysimum mediohispanicum, E. nevadense, and E. baeticum were collected from several populations located in the Baetic Mountains, South of Spain (Table 1 shows the code and location of all populations). Leaves were dried and preserved in silica gel until DNA extraction. Pre-opening flower buds at the same development stage were stored in liquid nitrogen for RNA extraction.
DNA extraction and sequencing. We used an individual sample for each species (Table 1) Chloroplast genome assembly and annotation. We assembled de novo the chloroplast genomes of E. mediohispanicum, E. nevadense, and E. baeticum using the NOVOPlasty pipeline v.6.2.3 7 (Fig. 1). Basically, through this pipeline a cp genome is assembled from whole genome sequencing (WGS) data, starting from a related single seed sequence iteratively extended bidirectionally until the circular genome is obtained. We used untrimmed reads as recommended by Dierckxsens et al. 7 and Arabidopsis thaliana cpDNA sequence (NC_000932.1) as the seed, since Erysimum is a close relative of Arabidopsis 39 . We specified the following parameters: automatic insert size detection, a genome range from 120000 to 200000, a K-mer value of 39, an insert range of 1.6, a strict insert range of 1.2, and the paired-end reads option. After assembling the full chloroplast genome of E. medihospanicum, we proceeded to assemble the cp genomes from the RNA-Seq data by using this chloroplast genome as a reference. From the RNA-Seq libraries, we first trimmed the adapters in the raw reads using cutadapt v.1.15 40 . For trimming adapters in 5′ and 3′ direction we used the "-b" option, and only used the prefix of the adapter sequence that is common to all "TruSeq Indexed Adapter" sequences (AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC). In addition, we used the "-n" option to search repeatedly for the adapter sequences (28 iterations). This option ensures that the correct adapters are detected by searching in loops until any adapter match is found or until the specified number of rounds is reached. Then, we quality-filtered the reads using Sickle v.1.33 41 , a trimming software which uses sliding-window analyses along with quality and length thresholds to cut and discard the reads which do not fit the selected threshold values. We specified the "pe" option for paired-end reads and the "-t" to use Illumina quality values (see https://github.com/najoshi/sickle). After filtering, we used the read mapper of Geneious R. 11 42 with the highest sensitivity and default parameters (http://www.geneious.com) 42 for a reference-guided assembly of the trimmed reads using the E. mediohispanicum reference assembly (see above). We validated the results obtained with Geneious by comparing them to read maps obtained with the BWA read mapper 43 .
The program cpGAVAS (Chloroplast Genome Annotation, Visualization, Analysis and GenBank Submission Tool 44 ) was used for annotating and visualizing the cp genomes. This program takes as input a FASTA file containing the genome information and performs bioinformatic analyses to annotate the genome. We used the annotated Arabidopsis thaliana cp genome (NC_000932.1) 45 . Protein coding genes were manually curated. Lastly, cpGAVAS gives as output the statistics of the annotation process, the annotated genome, and a visualization of the annotated genome. The annotations were then manually curated using Geneious R.11 42 . All transfer RNA sequences (tRNA) encoded in the cp genomes were verified using tRNAscan-SE 2.0 46 with the default search settings. The step-by-step process is presented in Fig. 1.

Comparative analysis among cp genomes assemblies.
To compare the cp genomes assembled from DNA or RNA libraries, we used the mVISTA software, part of the VISTA suite of tools for comparative genomics (http://genome.lbl.gov/vista/mvista/submit.html). This software compares DNA sequences from different species by pairwise alignment and allows the visualization of these alignments with annotation information. The output allows the identification of homologies between sequences, determining the percentage of identity between them using a sliding window of predefined length. We selected default parameters, a RankVISTA probability threshold of 0.5, and the Shuffle-LAGAN mode, which is a global alignment algorithm for finding rearrangements (inversions, transpositions, and some duplications). We used the A. thaliana cpDNA as a reference (NC_000932.1) 45 .
The sequence conservation profiles were visualized in mVISTA plots 47 .
We investigated the degree of within-genome variation of the assembled cp genomes. In particular, we performed a reference-guided assembly in which we remapped the quality-trimmed reads (as for the RNA-Seq assemblies, see above) to each assembled genome using the Geneious R. 11 42 mapper with medium-low sensitivity and default parameters (http://www.geneious.com) 42 . Later, we estimated the percentage of pairwise identity of each assembly. This statistic gives the average identity (as %), computed by scoring a hit when all pairs of bases are identical and dividing it by the total numbers of pairs.
For each species, we explored the degree of overall sequence variation found within the three replicas of RNA-Seq assembled genomes and then compared the results to those of a similar analyses that included also the genome assembled from genomic libraries. For this purpose, we estimated the nucleotide diversity (π) among the three replicas of cp genomes assembled from RNA-Seq, and then computed it again including the corresponding genomic library. Genomes were first aligned using MAFFT with the following parameters: FFT-NS-2 fast progressive method algorithm, a scoring matrix of 200PAM/k = 2, gap open penalty of 1.53 and offset value of 0.123. Then, we estimated the cpDNA nucleotide diversity using VariScan v.2.0.3 48 .
We studied the degree of sequence variation of some relevant chloroplast genes within the three replicas of RNA-Seq, and then explored it but including the genes assembled from genomic libraries. We first extracted and assembled all the chloroplast genes using the HybPiper pipeline v.1.2 49 . This pipeline uses BWA 43 to align reads to target sequences, and SPAdes 50 to assemble these reads into contigs. Once cpDNA genes were obtained, we selected 12 genes out of the total: rbcl, psaA, psbA, ndhK, atpA, atpH (with an important function in the photosynthesis process 51 ), rpoA, rps3, rrn16S, trnH (as self replication genes 52 ), yfc2 (the largest plastid gene in angiosperms 53 ), and matK (the only maturase of higher plants and widely used in angiosperm systematic 54 ). Then, we aligned these genes using MAFFT, as explained above. Lastly, we calculated the percentage of pairwise identity between the genes obtained from the three RNA-Seq replicas, and the same but including those from genomic libraries.
The size and location of repeat sequences, including palindromic, reverse and direct repeats, within these cp genomes were identified using REPuter software 55 . Following Asaf et al. 8 and Ni et al. 56 REPuter was parametrized with the following settings: Hamming distance of 3; 90% or greater sequence identity; and minimum repeat size of 30 bp.
Analysis of minimum transcriptome depth to produce quality cp genomes assemblies. To analyze the impact that sequencing depth has in the assemblage of transcriptome data into a complete cp genome, we subsampled the transcriptome reads of E. nevadense four times at 1 M, 5 M, 10 M, 20 M and 30 M paired reads. These reads were processed and mapped to the cpDNA of E. mediohispanicum with Geneious R.11 42 with medium-low sensitivity and default parameters as previously done. We calculated several mapping quality indexes (coverage of bases, expected errors, mean confidence, and % of Q40 positions) with Geneious R.11 42 and plotted them against the sub-sampling depth.
Cross-validation of the methodology. In order to estimate the recovery of complete cpDNA chromosomes from RNA-Seq libraries in other plant species, we downloaded five transcriptomes from the Sequence Read Archive website and processed them with our workflow. We downloaded two A. thaliana (SRR6757372; SRR6676021), one E. cheiri (SRR5195368), one Moricandia suffruticosa (SRR4296233), one M. arvensis (SRR4296231), one Oriza sativa (SRR7079258), and one Zea mays (ERR1407273) transcriptome. These libraries were trimmed and quality filtered using cutadapt v.1.15 40 and Sickle v.1.33 41 with the same parameters described above, and mapped using Geneious R.11 42

Results
Chloroplast genome assembly and annotation. Genomic DNA libraries. We assembled de novo the whole chloroplast genomes of three Erysimum species. The assembled genomes were circular and have a total length of 154,599 bp, 154,660 bp, and 154,581 bp in E. mediohispanicum, E. nevadense, and E. baeticum, respectively (Figs 2, S1 and S2). These chloroplast genomes displayed the typical quadripartite structure of most angiosperms (See Table 2 (Table S2). Thus, the number of unique protein-coding genes was 124 for the three species. These chloroplast genomes contained 29 unique transfer RNA genes and 8 unique ribosomal RNA genes (Fig. 3). The number of intra-gene regions was 150 in each cp genome. We found eight split genes (rpl2, atpF, rpoC1, psaA, ycf3, clpP, ndhB, ndhA; see Table 3) with intronic regions for each cp genome. Lengths of intronic regions are shown in Table S3. The overall GC content was 36.6%, indicating similar conserved GC levels among the Erysimum chloroplast genomes. A summary of the number of sequences assembled, mean assembly coverage, and percentages of pairwise identity are shown in Table S4 (Supporting Information).

RNA-Seq libraries.
We assembled the chloroplast genomes from three different replicas for each of the three species. We recovered high-quality complete chloroplast genomes that were very similar to those obtained from genomic DNA.  Table 2), as well as the protein-coding gene contents, tRNAs, and rRNAs were very similar between species replicates but when comparing with the chloroplast genomes obtained from genomic DNA, we found that the IRb regions were shorter and SSC regions were slightly larger (see Fig. S3 for chloroplast borders comparison). We found the same eight split genes with introns regions that were found in cp genomes obtained from genomic DNA (rpl2, atpF, rpoC1, psaA, ycf3, clpP, ndhB, ndhA; see Table 3). The lengths of all intronic regions are shown in Table S3. We found the same number of intra-gene regions using RNA-seq and genomic libraries (150 in each cp genome). The overall GC was 36%. A summary of assembly statistics including the bases assembled, mean assembly coverage, and percentages of pairwise identity are shown in Table S4. The results of the mapping assembly using BWA were highly similar (Table S5).
Repeat and SSRs analyses. The total number of repeats was 64, 78, and 65 in E. mediohispanicum, E.
nevadense, and E. baeticum, respectively. Forward repeats were the most common across the three species, followed by palindromic repeats. Reverse and complement repeats were found in low abundance (Table 3). In particular, E. mediohispanicum contained 38 forward, 24 palindromic, and two complement repeats; E. nevadense contained 51 forward, 25 palindromic, and one complement repeats; and E. baeticum contained 38 forward, 25 palindromic, and two complement repeats, respectively. In addition, the repeats from the three species had a sequence identity greater than 90%. The length of these repeats ranged for all the species from 30 to 26,429 bp, and the most common copy length had 30 bp. The number of repeats in the chloroplast genome assembled from RNA-Seq data was similar to that obtained from genomic DNA ( Table 3). The average number of repeats was  Table 2. Characteristics of the chloroplast genomes of Erysimum: type of library (genomic DNA or RNA-Seq library), length of the cp genome (bp), number of assembled reads, length of the two inverted repeats (IRa and the IRb), length of the small single copy (SSC), and of the large single copy (LSC) region, and GC% content. 78.3, 90, and 84, for the three replicas of E. mediohispanicum, E. nevadense, and E. baeticum, respectively. Forward repeats were the most common, followed by palindrome repeats, with lower levels of reverse and complemented repeats. The repeats of these population samples had a sequence identity greater than 90% for each species. The length of these repeats reached from 30 to 14, 353 bp, with the units with 30 bp being also the most common. The SSRs contained in the three chloroplast genomes were analyzed using the MISA Perl script (Fig. 4) Table S6 shows the numbers of SSR's that were quantitatively different between cp genomes assembled from genomic and RNA-Seq libraries. Among them, most of the SSRs were also mononucleotide repeats, with A/T repeats showing the highest proportion in the three replicas per species.
Genomic comparison. Results from mVISTA plots revealed a high similarity, with 99% of shared sequence identity in pairwise comparisons, between chloroplast genomes from genomic libraries and those from RNA-Seq libraries (See Fig. 5; the top and bottom percentage bounds are shown to the right of every row). These plots also showed a high degree of synteny between the three Erysimum species. In addition, the two IR regions were more similar than the LSC and SSC regions in all these species. Lastly, non-coding regions reveal a higher divergence than coding regions. Nucleotide diversity (π) was lower among the three replicas assembled from RNA-Seq. In contrast, nucleotide diversity increased dramatically (~ three orders of magnitude) when including the cp genome from genomic libraries in the alignments (0.35988 vs. 0.00008 for E. mediohispanicum; 0.36617 vs. 0.00037 for E. nevadense; and 0.36068 vs. 0.00123 for E. baeticum). Percentages of pairwise identity were always higher than 99% when comparing genes assembled from the different RNA-Seq replicas, and this similarity did not decrease when including genes assembled from genomic libraries (see Table S7).
Effect of sequencing depth. We assembled the chloroplast genomes from four different resampling of an E. nevadense transcriptome at 1 M, 5 M, 10 M, 20 M and 30 M paired reads. With this particular transcriptome, chloroplast genomes were obtained with coverage >95% from libraries of only 1 M reads, and with coverage >99% from sequencing depths >5 M reads (see supplementary information Figure S4). As expected, all metrics of mapping quality as well as the mean coverage at each position (p < 0.0001; R 2 = 0.921; from ~170 K to close to 200 K) increased significantly with sequencing depth (see supplementary information Figure S4).

Cross-validation results.
We assembled the cp genomes from RNA-Seq data of five species: A. thaliana, E.

Discussion
Our results showed that complete chloroplast genomes can be reliably assembled from transcriptomic data. We studied some Erysimum species as a proof-of-concept, and obtained genomes congruent in structure and sequence with previously published chloroplast genomes 2,3,58 . Both the chloroplast genomes assembled from transcriptomic and genomic libraries exhibit the typical quadripartite structure, low GC content, and are mainly composed of polythymine (polyT), and polyadenine (polyA) repeats 59 . Chloroplast genomes assembled from RNA-Seq data are highly similar in terms of SSRs, number of repeats, and plastome composition (CDS, exons, genes, rRNA, and tRNA) to those assembled from genomic libraries. Moreover, the similarity of the genomic and RNA-Seq assemblages validates that chloroplast genomes are fully transcribed. This is in line with findings from Shi et al. 34 who showed full transcription of chloroplast genome in photosynthetic eukaryotes using several tissues (flowers, complete seedlings, and seedlings shoots). Here, we show that chloroplast genomes of flower buds, the tissues we have used to obtain the RNA-Seq libraries, are also fully transcribed. Therefore, chloroplasts appear to be fully transcribed across organs and development stages in angiosperms, at least in samples containing functional plastids.
We have found significant differences in nucleotide diversity when comparing both kinds of assemblages (RNA-Seq vs genomic libraries). This may be explained by post-transcriptional modifications, i.e., by RNA-editing 60 . However, we found that nucleotide diversity greatly increased when including the assemblies from genomic libraries into the alignments. Accordingly, nucleotide diversity was lower when only comparing the three replicates of the RNA-Seq data. This implies that genomic assemblies were more heterogeneous or noisier than transcriptomic ones. Since both libraries were obtained using similar Illumina platforms, it appears that the genomic libraries were intrinsically more heterogeneous. This heterogeneity is likely caused by segments of chloroplast DNA transferred to the nuclear genome (i.e., nuclear plastic DNA or NUPT) that may potentially be incorporated during the mapping procedure introducing heterogeneity (i.e., within-genome polymorphism) into the cpDNA genomic assemblies 45,61 . However, NUPTs are generally fragmented and eliminated from the nuclear genome and therefore not transcribed, or transcribed at low level [62][63][64] , and therefore they should not be recovered in the RNA-Seq libraries. Moreover, the lack of differences in pairwise identity when only comparing genes from RNA-Seq to those from genomic libraries may be consequence of NUPTs located at the intergenic regions, as have been found in previous studies (e.g., only 25% of NUPTs in Arabidopsis thaliana are located in genes 65 ). NUPTs are well documented in plants 66 , and they often represent a significant part of the nuclear genome 67,68 . Because of the maternal inheritance in most plant genera 69 , cpDNA is widely used for the inference of relationship among plants. Therefore, the presence of NUPT into cp genomes may lead to erroneous phylogenetic inferences 66 . According to our results, using cpDNA assembled from transcriptomes might reduce the problems due to NUPT inclusion when using cpDNA in phylogenomics. Alternatively, methods specifically designed to correct these assembly errors have been developed for genomic data, such as the dnaLCW method 61 , and should be consider whenever possible. However, validating that NUPTs are a source of error in cp genome assembly requires comparison with a reference genome, which is currently not available for Erysimum. Therefore, the potential misleading mapping caused by this type of genetic elements will require further studies.
When we tested our methodology across several plant species, we found that assembling chloroplast genomes from RNA-Seq data is a relatively fast and flexible approach. In the light of these results, we put forward a pipeline-like procedure in the hope that it can be useful to other researchers (Fig. 1). In addition, we showed that, although the chloroplast genome coverage increased with the number of reads used for the assembly, 1 M reads was sufficient to obtain a 95% coverage of the cp genome. These results corroborate that the chloroplast could be fully transcribed and is easily assembled from transcriptomic data even at low-medium coverage. Moreover, cross-validation (Supplementary Table S8) showed that assembling the cp genome using transcriptomes from the SRA database is feasible even though the reliability of the assemblage is always a function of the tissue and methodology used. For example, Arabidopsis thaliana cp genomes, assembled from RNA-Seq data coming from different libraries (SRR667021 and SRR6757372), produced different assembly results that were related to differences in coverage and number of reads. Furthermore, the genome of E. cheiri cpDNA was surprisingly not fully assembled despite being a closely related species to the Erysimum species used in this study 37 . However, this result may be explained by the fact that this E. cheiri transcriptome was obtained from petals. The reliability of our results is probably attributable to careful sample preparation (our RNA-Seq samples were submitted to a treatment that depleted rRNA implying that the samples were enriched in the other types of RNA), and because sequencing depth, at least over a minimum threshold ~5 M reads (see Figure S4), does not appear to be a crucial factor. Therefore, as a general rule, samples obtained from photosynthetic tissues, depleted in rRNA and high quality sequenced (as indicated by quality scores) are likely to be trustworthy. We conclude that assembling cp genomes from good quality transcriptomic data (either obtained de novo or downloaded from public databases such as the SRA database) may be a straightforward approach in plant systematics and phylogeny. In fact, this approach may reduce the risk of incorporating NUPTs avoiding posterior phylogenetic incongruences, although precautions must be taken due to the possibility of RNA editing, and alternative methods 61 could be also used to minimize the assembly of NUPTs, or other nuclear DNA, into cp genomes. In summary, we think the pipeline presented here is an accessible and time saving approach to produce high-quality cp genomes that could complement other genomic approaches.