Revisiting the Zingiberales: using multiplexed exon capture to resolve ancient and recent phylogenetic splits in a charismatic plant lineage

Taxon sampling, DNA extraction, and library preparation

Sampling included several members of each of the eight families: Heliconiaceae (5), Musaceae (9, including 2 previously published whole genomes, D’Hont et al., 2012; Davey et al., 2013), Strelitziaceae (3), Lowiaceae (2), Zingiberaceae (16), Costaceae (10), Marantaceae (7), and Cannaceae (3). In total, 53 individuals were sequenced de novo (Table S1). DNA was extracted using an SDS and salt extraction protocol (Edwards, Johnstone & Thompson, 1991; Konieczny & Ausubel, 1993) from freshly collected leaves dried in silica, eluted in TE buffer, and sonicated with a Bioruptor® (Diagenode, Liège, Belgium) or qSonica Q800R machine to an average size of approximately 250bp. Sonicated DNA was cleaned and concentrated with solid phase reversible immobilization magnetic beads (Sera-Mag; GE Healthcare, Little Chalfont, UK), and libraries were prepared according to Meyer & Kircher (2010).

Probe design, sequence capture, sequencing

To generate a nuclear probe set, the Musa acuminata CDS was downloaded from the banana genome hub (http://banana-genome.cirad.fr/) and split into annotated exons. Raw reads of transcriptomes for each of the remaining seven families were cleaned to remove adapters, low-complexity sequences, contamination, and PCR duplicates (Singhal, 2013). Cleaned transcriptome reads were aligned to the Musa acuminata exons using NovoAlign v3.01 (http://novocraft.com) with –t 502 to allow highly divergent sequences to map. After mapping, SNPs were called using SAMtools v0.1.18 (Li et al., 2009) and VarScan v2.3.6 (Koboldt et al., 2012) and consensus sequences for each family were made based on SNP calls. All exons were filtered for: (1) having overlapping read coverage in all 7 families (2) being longer than 150 bp (3) having between 30–70% GC content (4) being unique by reciprocal BLAST (5) not being found in the RepeatMasker database (command parameters can be found in Supplementary Methods). Although conservation to the Musa sequence was not used as a direct filter, the alignment protocol inherently limits chosen regions to those with relative conservation across the families–the minimum percent identity between any Musa exon and a family specific bait for the same region was 86% (Table S1). After filtering, a total of 494 exons from 418 genes for each of the eight families (the Musa reference sequence plus each sequence from the seven families) were printed with 1 bp tiling twice each on an Agilent 1M microarray chip (G3358A) (Fig. 1A). A second chip was printed with one complete plastid genome from each family (Barrett et al., 2014) with slightly less than 1 bp tiling. Libraries from a total of 56 individuals were quantified by Qubit® and pooled in equimolar quantities. The total library pool was split in half and one half was hybridized to the nuclear array and the other half was hybridized to the plastid array (Hodges et al., 2009). After hybridization, pools were subject to a limited amount of PCR amplification and enrichment success was verified with qPCR using primers matching both targeted and non-targeted regions. Because of known bias toward plastid dominance in sequenced reads owing to a greater percentage of plastid DNA in the total genomic DNA extractions, the separate hybridization pools were combined in a ratio of 3 parts nuclear to 1 part plastid and sequenced (100 bp paired-end reads) in one lane of a Illumina® HiSeq® 2500 platform at the Vincent J. Coates Genomics Sequencing Facility at the University of California, Berkeley.

Read processing

Raw reads were cleaned to remove adapters, low-complexity sequences, contamination, and PCR duplicates (Singhal, 2013). Custom Perl scripts were created to perform a series of alignment and reference adjustments using NovoAlign v3.01 (NovoCraft: http://novocraft.com), VarScan v2.3.6 (Koboldt et al., 2012) and Mapsembler2 v2.1.6 (Peterlongo & Chikhi, 2012) to generate a per individual reference for SNP calling without the need for de novo assembly (Fig. 1B). Perl scripts are available in a github repository (https://github.com/chodon/zingiberales). The plastid sequences were processed the same way except extension with Mapsembler2 was omitted, and individual genes were extracted from the whole plastid prior to final mapping. Finally, reads were mapped with NovoAlign with –t 90 and PCR duplicates were removed with Picard v1.103 (http://picard.sourceforge.net). SNPs were called following best practices guidelines using the HaplotypeCaller and readBackedPhasing algorithms in GATK v3.1.1 (McKenna et al., 2010; DePristo et al., 2011; Van der Auwera et al., 2013), except quality scores were not recalibrated because the lack of a reference set of known variants. Consensus sequences were created based on SNP calls for regions with greater than 20× coverage (Nielsen et al., 2011). SNPs in areas with less than 20× coverage were converted to Ns and regions with less than 5× coverage were discarded. For outgroup taxa, raw reads from transcriptomes generated as part of OneKP were subject to the same pipeline as sequences generated de novo. The raw sequence data from the Musa balbisiana genome project (Davey et al., 2013) was also subject to the pipeline, but only aligned for the plastid gene set. Raw de novo sequence reads are available on NCBI under bioproject SRP066318 and the final concatenated alignment is accessible on github https://github.com/chodon/zingiberales.

Alignment

After consensus sequences were made, a second pipeline was made to pass sequences through a series of alignment steps to (1) trim sequences to the Musa reference (MAFFT v7.164 (Katoh et al., 2002; Katoh & Standley, 2013) and mothur v1.34.4 (Schloss et al., 2009)), (2) place sequences into coding frame (MACSE v1.01b (Ranwez et al., 2011)), and (3) align by codon position (prank v140603 (Löytynoja & Goldman, 2005)). Plastid gene introns were spliced out by hand in Geneious v5.6.4 (Kearse et al., 2012) prior to step 3, above. After alignment, several additional steps were taken to eliminate genes that might contain non-orthologous sequences. Gene trees were generated with RAxML v8.1.17 (Stamatakis, 2014) and the single gene trees were assessed to identify those in which the gene of a single individual taxon accounted for greater than 15% of the total tree length (dos Reis et al., 2012). Exon sequences from one individual were BLASTed to the nucleotide collection database (BLASTN v2.2.30+, Altschul et al., 1997). Exons were removed from further analyses if significant BLAST hits were found to a whole plastid genome, or to ribosomal, transposon, or mitochondrial DNA. Exons were also removed from further analysis if they had unexpectedly high average coverage of greater than 200× or because frameshifts were introduced during codon position assignment or the alignment had too many indels to be reliable (Table S2). We also manually checked all alignments for potential problems (Rothfels et al., 2015). Command parameters for all steps can be found in Supplementary Methods.

Phylogenetic analyses

The nuclear and plastid sequence data were concatenated and analyzed using maximum parsimony (MP), maximum likelihood (ML), and coalescent approaches. For MP, PAUP* v4.0a142 (Swofford, 2002) was used to perform a heuristic search with 100 random addition sequence replicates and default parameters (TBR branch swapping with one tree held per replicate). MP support was evaluated with 1000 bootstrap replicates, each with 10 random addition sequence replicates. For ML reconstruction, gene-by-codon position partitions were created for the complete concatenated data set resulting in a total of 1128 initial partition subsets. These initial subsets were then grouped using the relaxed hierarchical clustering algorithm with a 1% search strategy (Lanfear et al., 2014) implemented in PartitionFinder v1.1.1 (Lanfear et al., 2012). The resulting partitioning scheme generated by PartitionFinder consisted of 112 subsets (see Supplementary Methods). The PartitionFinder scheme was analyzed with RAxML v8.1.24 (Stamatakis, 2014) with the GTR+Γ₄ model of sequence evolution estimated for each partition subset and the topology linked across partitions. ML support was evaluated for the same partitioning scheme with 1000 bootstrap replicates, using the rapid bootstrap algorithm (Stamatakis, Hoover & Rougemont, 2008), and using the CAT₂₅ approximation instead of Γ₄, to model site-to-site rate heterogeneity (Stamatakis, 2006). The RAxML analysis was performed on the CIPRES web server (Miller, Pfeiffer & Schwartz, 2010). For coalescent analysis, first, best ML gene trees and bootstrap gene trees were generated in RAxML v8.1.17 using the GTR+Γ₄ model. Best gene trees were the best tree of 20 independent searches and support was evaluated with 1000 bootstrap replicates. ASTRAL-II v4.7.8 (Mirarab & Warnow, 2015) was used to generate an optimal tree based on the best ML gene trees and support evaluated with the 1000 bootstrap replicates.

Probe design, sequence capture, and alignment

All targeted regions for all individuals were successfully captured, although average coverage varied based on gene region (Fig. 2A), individual, and phylogenetic distance to the reference sequence (Fig. 2B). Members of the Musaceae, in general, captured better than any other family, likely because they are phylogenetically closest to the original genomic reference upon which the probes were designed. Within each family, close relatives of the species or taxon used to design the bait had higher success rates of capture than more distant members of the family. For example, Siphonochilus kirkii, had the lowest average coverage and capture efficiency for Zingiberaceae (Figs. 2B and 2C) as predicted by its evolutionary distance from the transcriptome-sequenced taxon Curcuma longa. The minimum percent identity of any captured sequence to its bait was 73%, while average distances were between 94–99% identity (Table S1). Of the total sequenced bases, the capture efficiency varied across individuals with the maximum percentage of bases mapping 3.5× higher than the minimum percentage (Fig. 2C). An average of 26% of captured bases mapped to target, which is similar to capture efficiency reported in captures of human mitochondrial DNA (Maricic, Whitten & Pääbo, 2010) and transcriptome based capture of chipmunk DNA (Bi et al., 2013). Despite the attempt to capture nuclear and plastid targets evenly, sequencing was highly biased towards plastid targets (Fig. 2C). There was some variability between individuals that was independent of phylogenetic distance, likely due to the standard variation in the success of DNA library preparation, which results from differences in DNA quality, genome size, and difficulties of accurately quantifying DNA for pooling in equimolar quantities. Any differences in DNA concentration were likely amplified in the post-hybridization PCR enrichment step.

Of the 494 nuclear probe exons, 124 were removed from further analyses based on coverage, BLAST results, skewed tree length, or alignment anomalies (Table S2). These 124 exons were from 110 genes. Twenty exons from 14 genes failed a test for coverage outliers because they had greater than 200× coverage, which is outside of the 99.99% confidence interval (Fig. S1). It is possible that these regions were either incorrectly annotated as nuclear regions in the Musa draft genome, or were transferred to the nuclear genome from more high copy genomes, especially considering that 15 of these exons were annotated as having an “unknown chromosomal location” in the Musa draft genome (Fig. 2A, Table S2). A total of 37 exons from 34 genes were removed from the nuclear dataset and 13 genes from the plastid dataset due to skewed tree length. Four nuclear exons from two genes were removed because of introduced frameshifts and ycf1 from the plastid was eliminated because of insertions and deletions in the alignment apparent after manual inspection. Finally, 63 additional exons from 61 genes were removed because of a top BLAST hit to a whole plastid genome, mitochondrial, transposon or ribosomal DNA. Of these 63 exons, the 27 ribosomal and 21 mitochondrial exons could likely be included in further analyses or within family specific analyses in future work after analyzing secondary structure and genomic location.

The final dataset of 308 nuclear genes had a total aligned length of 81,546 bp with 24,379 (29.9%) parsimony informative sites and an average coverage of 40 ± 13× (mean ± s.d.). The 68 gene plastid dataset had a total aligned length of 56,202 bp with 8,336 (14.8%) parsimony informative sites and an average coverage of 377 ± 589× (Table S2).

Phylogenetic analyses

The recovered topology (Fig. 3) places Musaceae as sister to all other families with 100% parsimony bootstrap support (pb) and maximum likelihood bootstrap support (mlb). The ginger families (Cannaceae, Costaceae, Zingiberaceae and Marantaceae) are well supported (100 pb/mlb) as monophyletic. The MP and ML trees are largely congruent and support values are generally high from shallow to deep phylogenetic relationships. The optimal coalescent tree does not conflict with the ML tree but there is low support for several nodes, notably the placement of Musaceae sister to the rest of the order (Fig. 3).