Fresh leaf samples of six Amana species, five from China and one from Japan (Table 1), were field-collected and dried with silica gel. Voucher herbarium specimens were deposited at the Herbarium of Zhejiang University (HZU). We extracted total DNA from ca. 3 mg of the silica-gel dried leaf tissue for each species using DNA Plantzol Reagent (Invitrogen) and following the manufacturer's protocol. The qualities and quantities of genomic DNA were checked on an Agilent BioAnalyzer 2100 (Agilent Technologies). Short-insert (500 bp) paired-end libraries were generated by using Genomic DNA Sample Prep Kit (Illumina) according to the manufacturer's protocol. Genomic DNA of each species was indexed by barcode tags and then pooled together for sequencing in one lane of HiSeq™ 2500 (Illumina, San Diego, California, USA) at Beijing Genomics Institute (BGI, Shenzhen, China).

For each Amana species, raw reads (125 bp read length) were firstly cleaned by removing low-quality reads with Phred scores of <20 using the CLC-quality trim tool (http://www.clcbio.com/products/clc-assembly-cell/). Secondly, we assembled the clean reads into contigs on the CLC de novo assembler (http://www.clcbio.com/products/clc-assembly-cell/), under the following settings: minimum contig length of 200 bp, mismatch cost of 2, deletion and insertion costs of 3, length fraction of 0.8, and similarity fraction of 0.8. Thirdly, all the contigs were aligned to the reference genome (Lilium longiflorum Thunb., KC968977) using BLAST (http://blast.ncbi.nlm.nih.gov/), and aligned contigs were oriented according to the reference genome. Then, contigs were aligned with the reference genome for constructing the draft chloroplast genome of each Amana species in GENEIOUS V9.0.5 (http://www.geneious.com). Finally, clean reads were re-mapped to the draft genome and yielded the complete chloroplast genome sequences.

Preliminary annotation of these Amana chloroplast genomes was conducted on the program Dual Organellar GenoMe Annotator (DOGMA; Wyman et al., 2004). DOGMA annotations were further corrected for the start/stop codons and intron/exon boundaries by comparison with homologous genes from L. longiflorum (GenBank accession no. KC968977) and Fritillaria hupehensis (GenBank accession no. KF712486) using MAFFT v7 (Katoh and Standley, 2013). In addition, tRNAscan-SE (Schattner et al., 2005) was used to verify the tRNA genes with default parameters, and the ultimately annotated chloroplast genomes were deposited in GenBank (accession numbers listed in Table 1). The cp genome maps were drawn in OrganellarGenome DRAW (Lohse et al., 2007). Codon usage, as well as relative synonymous codon usage (RSCU, Sharp and Li, 1987) value was estimated for all exons of protein-coding genes with the program CODONW V1.4.2 (http://codonw.sourceforge.net/).

Chloroplast genome comparisons across the six Amana species was performed in Shuffle-LAGAN mode on the mVISTA program (genome.lbl.gov/vista/index. shtml, Frazer et al., 2004), using the annotation of A. kuocangshanica as a reference. To evaluate whether different cp genome regions underwent different evolution patterns in this genus, and to explore highly variable regions for future population genetic and species identification studies, we sequentially extracted both coding regions and noncoding regions (including intergenic spacers and introns) after alignment using MAFFT v7 under the two criteria that aligned length is >200 bp and at least one mutation site is present. After that, the nucleotide variability of these regions was evaluated with DNASP V5.10 (Librado and Rozas, 2009).

REPUTER (Kurtz and Schleiermacher, 1999) was used to determine the size and position of repeat sequences, which included direct, inverted, complement and reverse repeats in the Amana chloroplast genomes. The minimum length of repeat size and sequence identity was set to 30 bp and >90%. MISA perl script (Thiel et al., 2003) was applied to detect the simple sequence repeats (SSRs) in the six Amana cp genomes with thresholds of 10, 5, 4, 3, 3, 3 repeat units for mono-, di-, tri-, tetra-, penta-, and hexanucleotide SSRs, respectively.

Phylogenetic analyses were conducted on the six Amana species and one species each for Erythronium, Tulipa and Lloydia, using L. longiflorum (KC968977) and Fritillaria cirrhosa D. Don (KF769143) as outgroups based on previous studies (Rønsted et al., 2005; Kim et al., 2013). Chloroplast sequences of these 11 species were aligned using MAFFT v7. In order to evaluate possible alternative hypotheses of phylogeny, topologies were constructed by both maximum likelihood (ML) and Bayesian inference (BI) methods using not only the complete cp genome sequences (162,505 bp), but also the exons of protein-coding genes (78,815 bp). We also tried two different partitioning strategies for the second dataset: (1) separating each gene as a partition, (2) divided the data matrix into three partitions, corresponding to the first, second and third codon positions.

The best-fitting models of nucleotide substitutions were determined by the Akaike Information Criterion (AIC) in JMODELTEST V2.1.4 (Posada, 2008). The GTR+I+G model was most suitable for both datasets. Maximum likelihood analyses were conducted using RAXML-HPC v8.2.8 (Stamatakis, 2014) with 1000 bootstrap replicates at the CIPRES Science Gateway website (Miller et al., 2010). Bayesian inference (BI) analyses were performed in MRBAYES v3.2 (Ronquist and Huelsenbeck, 2003). Two independent Markov Chain Monte Carlo chains were calculated simultaneously for five million generations with trees sampled every 500 generations. The first 25% of calculated trees were discarded as burn-in, and a consensus tree was constructed from the remaining trees to estimate posterior probabilities (PPs).

After filtering the low-quality reads and adaptor sequences, 12,248,447–28,408,770 clean reads (of 125 bp length) were obtained for the six Amana species. Through de novo assembly, 10,585 contigs (A. latifolia) to 70,038 contigs (A. anhuiensis) were assembled with N50 contigs varing from 338 to 348 bp (Table 1). Subsequently, three to six initial contigs which were found to be significantly homologous to the reference genome were combined to generate each chloroplast genome, with no gaps or missing nucleotides (Ns) found.

The full length of the six Amana chloroplast genomes ranged from 150,613 to 151,136 bp (Table 1; Figure 2). Akin to other angiosperms, chloroplast genomes of the six Amana species present a typical quadripartite structure, including a pair of inverted repeat regions (IR with 25,629–25,859 bp) separated by one large single-copy region (LSC with 81,482–82,218 bp) and one small single-copy region (SSC with 17,366–17,465 bp). All the complete cp genomes with annotation were deposited in GenBank (accession numbers listed in Table 1). The GC content in the LSC, SSC, and IR regions, and also in the whole genome sequences, were nearly identical among the six Amana species (Table 1).

The six Amana chloroplast genomes contained the same 132 genes, of which 20 were duplicated in the IR regions and 112 were unique genes comprising four rRNA genes, 30 tRNA genes and 78 protein-coding genes (Table 2). Of the 112 distinct genes, 15 held a single intron (nine protein-coding genes and six tRNA genes) and three (ycf3, clpP, and rps12) possessed two introns. The gene infA was lost in all six Amana species. The gene rps12 was trans-spliced; the exon at the 5′ end was located in the LSC region, however the 3′ exon and intron were located in the IR regions. The regions ycf15, ycf68, and ycf1 were identified as pseudogenes because they contained several internal stop codons. Besides, the rps19 gene located in the IRa/LSC junction region lost their protein-coding ability because of incomplete gene duplication. The similar event was also observed in the ycf1 region at the IRb and SSC border (see below for detailed information).

The number of codons encoded by 78 protein-coding genes in the six Amana chloroplast genomes ranged from 24,349 to 24,380. Amana kuocangshanica was randomly selected as an example for detailed investigation, owing to the similar result of codon usage and RSCU values for these six species (Table S1). Among the 24,369 codons in A. kuocangshanica, 2,529 (10.38%) encoded leucine and 286 (1.17%) encode cysteine, which were the most and least frequent amino acids, respectively (Table S1, Figure S1). In synonymous codons, RSCU value increased with the number of the codons. Furthermore, the RSCU values of 30 codons were greater than 1, suggesting that they are biased codons in the A. kuocangshanica chloroplast protein-coding genes (Table S1, Figure S1).

The IR/SC borders with full annotations for adjacent genes were compared among the six Amana chloroplast genomes (Figure 3). Except for A. anhuiensis and A. erythronioides, all the IRb regions expanded by 104–106 bp toward the rps19 gene with corresponding pseudogene fragment ψrps19 created at the IRa/LSC border. Long ψycf1 fragment with 1,121–1,154 bp was located at the IRb regions because the border between SSC and IRa extended into the ycf1 genes. In addition, the ndhF gene in A. wanzhensis overlapped with the IRa/SSC border by 66 bp. However, for the other five Amana species, the distance between ψycf1 and ndhF varied from 5 to 42 bp (Figure 3).

We analyzed the comprehensive sequence divergence of the six Amana cp genomes using the mVISTA software with the annotation of A. kuocangshanica as a reference. A genome-wide alignment revealed globally high sequence similarity (>90% identity) among them (Figure S2). Inverted repeat regions show a lower level of sequence divergence than LSC and SSC regions. In addition, 120 regions were eventually extracted to calculate the nucleotide variability, and the Pi value ranged from 0.02% (rrn23s) to 1.66% (rps15–ycf1). A total of nine regions (rps15–ycf1, accD–psaI, petA–psbJ, rpl32–trnL, atpH–atpI, petD–rpoA, trnS–trnG, psbM–trnD, and ycf4–cemA) with a nucleotide diversity >0.9% were recognized as hotspot regions that could be developed as molecular markers for future phylogenetic analysis and plant identification studies (Figure 4, Table S2).

A total of 371 repeats including 161 forward, 195 palindromic and 15 reverse repeats were identified in the six Amana cp genomes using the REPUTER software. Amana kuocangshanica possessed the greatest total number of repeats (69), while A. anhuiensis contained the fewest (54) (Figure S3A). For each Amana species, the majority of repeats (62.3% in A. kuocangshanica – 76.2% in A. edulis) ranged in size between 30 and 40 bp (Figure S3B). Repeats located in homologous regions with the same lengths were identified as shared repeats. Under this criterion, 38 repeats were shared by all Amana species and seven repeats were shared by five of the Amana species (all except A. edulis). Amana edulis showed the most distinct repeats (19), whereas A. anhuiensis had the least number (2) (Figure S3C, Table S3).

Each Amana species contained 69 (A. latifolia) to 76 (A. wanzhensis) SSRs, and more than half were composed of A or T bases (Figure S4A, Table S4). In the total 438 SSR regions, the proportion of the repeats situated in the intergenic spacer (IGS) regions led to 58.68%, while the regions located in the coding DNA sequence (CDS), CDS introns, tRNA introns and ψycf1 accounted for 15.53, 14.38, 4.34, and 7.08%, respectively (Figure S4B, Table S4). In addition, 29 SSRs (excluding mononucleotide SSRs) were identified as polymorphic SSRs between Amana species (Table S5), which could be useful for further population studies. Three criteria for identification were followed: (1) SSRs possessed the same repeat units (2) the number of repeat units is different and (3) SSRs located in the homologous regions.

In the present study, two datasets (whole chloroplast genome sequences and shared protein-coding genes) from the six Amana species, together with Erythronium sibiricum, Tulipa altaica, Lloydia tibetica and two outgroups, were used to conduct various phylogenetic analyses. Both ML and BI methods, based on different datasets and partitioning strategies, produced highly congruent topologies (Figure 5). The phylogenetic trees based on complete genome sequences, which had full support at every nodes [ML bootstrap (BS) = 100%, Bayesian posterior probabilities (PP) = 1], are shown here. Amana, Erythronium, and Tulipa were fully supported as a monophyletic group, in which the Amana species were resolved as a monophyletic clade that was sister to E. sibiricum. For the six Amana species, A. edulis was sister to a clade of the other remaining species with (A. erythronioides + A. kuocangshanica) and A. latifolia sharing a common ancestor and sister to (A. anhuiensis + A. wanzhensis) (Figure 5).

Our results revealed that the overall gene content and arrangement within the six Amana cp genomes are largely similar. This is expected considering the morphological similarities and presumed recent age of divergence among them. At the same time, however, we have been able to document that the plastomes of these species do, indeed, vary, even if their differences are small. The IR/LSC boundaries in the Amana chloroplast genomes (except those of A. erythronioides and A. anhuiensis) expand into the rps19 gene. This is congruent with a typical monocot cp genome structure (Yang et al., 2013). Among other taxa in the family Liliaceae, IR expansion into rps19 has been observed in Lilium (Kim and Kim, 2013), Fritillaria (Li et al., 2014), and Cardiocrinum (Lu et al., 2016). This suggests that the expansion of the IR/LSC junctions into rps19 may be an ancestral symplesiomorphy of the family Liliaceae, and thus provides no relevant phylogenetic information for addressing intrafamilial questions.

Millen et al. (2001) suggested that infA, which codes for translation initiation factor 1, has been entirely lost or has become a pseudogene approximately 24 separate times in 309 angiosperms. According to their results, the parallel loss of infA from the chloroplast genome occurred in both Tricyrtis (Liliaceae) and Smilax (Smilacaceae), which are members of Liliales. Kim and Kim (2013) revealed that this event also occurred in Alstroemeria (Alstroemeriaceae), which is closer to the basal Liliales than Smilax. Besides, they found that infA existed but seemed to have lost its function in the Lilium (Liliaceae) cp genome, because it had AAT instead of ATG in the start codon position and includes two premature stop codons. The pseudogenization of infA was also found in Fritillaria (Li et al., 2014) and Cardiocrinum (Lu et al., 2016), close relatives of Lilium. However, our study indicates the loss of infA in Amana. Overall, it shows that infA may have been lost in the most recent common ancestor (MRCA) of Liliales, and the pseudogenization of infA seems to be a synapomorphy of the (Lilium + Fritillaria) + Cardiocrinum clade. Further study is needed to improve our understanding of infA gene evolution in Liliaceae.

Phylogenetic analyses of the complete chloroplast genome sequences and separate analyses restricted only to 74 common plastid protein-coding genes both produced a well-resolved phylogenetic tree (Figure 5). The close relationship and relatively short branches among Amana, Erythronium, and Tulipa was confirmed, which is congruent with previous studies (Hayashi and Kawano, 2000; Allen et al., 2003; Rønsted et al., 2005; Zarrei et al., 2009; Clennett et al., 2012; Kim et al., 2013; Petersen et al., 2013). However, our phylogenetic trees unambiguously revealed a sister relationship between Amana and Erythronium (BS = 100%, PP = 1), clearly supporting the separation of Amana from Tulipa, as others have argued (Rønsted et al., 2005; Zarrei et al., 2009; Clennett et al., 2012; Kim et al., 2013). Although a closer relationship between Amana and Tulipa, or of Erythronium and Tulipa, has been suggested by some former studies, most of them used only a single locus and/or found only weak to medium support for their topologies (rbcL: Hayashi and Kawano, 2000; matK: Allen et al., 2003; ITS: Zarrei et al., 2009). In at least one case (Christenhusz et al., 2013), although a multi-locus phylogenetic analysis of these three genera showed Erythronium sister to Tulipa rather than to Amana, no external outgroup taxa from the sister clade were used to orient the tree, and so that the relative branching order of Amana, Erythronium, and Tulipa was actually undetermined in that study. The chloroplast genome gives the most strongly supported indication of relationships among the three genera, but the possible concordance between plastid gene trees and species trees remains tentative, given that the chloroplast genome sequences of ~150 kbp still essentially represent a single-locus (linkage group) phylogeny (Ruhsam et al., 2015).

Within Amana, the rare species may be recently evolved ecotypes of the widespread A. edulis, quickly adapting lineage. As ecotypes, these rare species are expected to show no genetic distinction at neutral loci, and may not merit species recognition lineage (Oberle and Schaal, 2011). Nevertheless, the six Amana species exhibit sequence divergences in plastid genomes (Figure S2), ruling out the possibility of being ecotypes. Furthermore, our phylogenomic analyses as expected recovered Amana as a monophyletic genus (BS = 100%, PP = 1), and strongly supported its division into two clades: a widespread species (A. edulis) and a clade of five rare species (BS = 100%, PP = 1; Figure 5). The result indicates that although they are allopatrically distributed across East China/South Japan, the five rare species share a more recent common ancestor with each other than they do with A. edulis. Therefore, it is unlikely that any of these rare species originated from the widespread species through local geographical and ecological isolation by progenitor–derivative speciation (Crawford, 2010). In fact, the two sister lineages exhibit different eco-geographies: while A. edulis is widespread in lowland evergreen broad-leaved or temperate deciduous forests of East/North China, Japan, and the Korean peninsula, the rare species are endemic to the montane warm-temperate-deciduous (WTD) forest in East China/South Japan. In line with evidence from palaeomodeling of East Asian forest biomes (Harrison et al., 2001) and recent phylogeographic studies (reviewed in Qiu et al., 2011), the exceptionally high diversification rate in the “rare-species” clade is mainly driven by long-term allopatric population isolation (viz. vicariance) in which climate-induced eco-geographic isolation through (a)biotic displacement of WTD forested habitats at different spatial–temporal scales and over glacial and interglacial periods is the primary vicariance factor (see also Qiu et al., 2009). Overall, our phylogenomic analyses based on chloroplast genomes have provided the first successful attempt to clarify intrageneric relationships within Amana. However, based on distributional considerations, hybridization is still expected to occur between the widespread A. edulis and rare species within their zone of sympatry. Although these cp genome data have generated a fully resolved phylogeny of the genus Amana (Figure 5), it is not possible to use such data to classify hybridization events because cpDNA is generally uniparentally inherited (Birky, 1995). In the future, multi-locus phylogenies, phylogeography and palaeo-climatic niche modeling are required to explore the time scales and demographies of species divergences as well as hybridization in this genus.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.