Plastome Genome Sequencing and Assembling

All genome data were sequenced using an Illumina Hiseq 2500 platform by Biomarker Technologies, Inc (Beijing, China). High-quality reads were obtained using the CLC Genomics Workbench v7.5 (CLC Bio, Aarhus, Denmark). Reference-guided assembly was then performed to reconstruct the chloroplast genomes using the program MITObim v1.7 (Christoph et al., 2013). In order to obtain accurate sequences, each species was assembled four times with the reference genomes A. cepa (KM088014), A. sativum (KY085913), A. victorialis (NC_037240), and A. obliquum (LT699701). Gaps that appeared in the assembled cp genomes were corrected by Sanger sequencing and the primers were designed using Lasergene 7.1 (DNASTAR, Madison, WI, United States). The primers and amplifications were shown in Supplementary Table S2. The program DOGMA (Wyman et al., 2004) was used to annotate the whole cp genome, and subsequently corrected within GENEIOUS R11 (Biomatters, Ltd., Auckland, New Zealand). Final plastid genome maps were drawn using OGDRAW (Lohse et al., 2013).

GC Content and Species Pairwise Ka/Ks Ratios

GC content of the complete chloroplast genome (CCG) and the third position GC content of codons in each species were calculated using PAML v4.8 (Yang, 2007). Each CDS sequence of all species was extracted and aligned with MAFFT v. 7 (Katoh and Standley, 2013). Pairwise Ka/Ks ratios of all species were calculated using the concatenated single-CDS alignments with KaKs Calculator version 2.0 (Wang et al., 2010).

SSRs Characterization and Chloroplast Genome Nucleotide Diversity

Perl script MISA (Thiel et al., 2003) was used to search microsatellites loci in the cp genomes with parameters being set as 10, 5, 4, 3, 3, and 3 for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotides, respectively. The DnaSP version 5.1 (Librado and Rozas, 2009) was used to calculate the nucleotide diversity of genes in LSC, SSC, and IR regions.

Indices of Codon Usage

Codon usage in these genes was assessed using the program codon W 1.4.4 (J. Peden)¹. Five values were used to estimate the extent of bias toward codons: the codon adaptation index (CAI), codon bias index (CBI), frequency of optimal codons (Fop), GC content of synonymous third codons positions (GC3s), and the effective number of codons (ENC).

Phylogenetic Analyses

In order to investigate the relationships of the six Allium species, all available complete genome sequences in allied families were downloaded from NCBI, including 20 species from Asparagaceae, 10 species in Allioideae (all Allium species), one species in Agapanthoideae, 2 species in Asphodeloideae, one species in Asphodelaceae and one species in Iridaceae (Supplementary Table S3). Firstly, all single-copy genes were extracted from all 41 taxa, and alignments of each gene were generated and trimmed. These alignments were then concatenated to produce an alignment of all single copy genes, which were used for phylogenetic analysis. Maximum likelihood (ML) analyses were performed using RAxML 8.2.8 (Stamatakis, 2014) with GTR + G model and 1,000 bootstrap replicates. Bayesian analyses were performed with MrBayes v. 3.2.5 (Ronquist and Huelsenbeck, 2003) under the GTR + I + Γ substitution model. The Markov chain Monte Carlo (MCMC) algorithm was run for 1 × 10⁸ generations, with one tree sampled every 1000 generations. The first 20% of trees were discarded as burn-in, and the remaining trees were used to build a 50% majority-rule consensus tree. The stationarity was considered to be reached when the average standard deviation of split frequencies remained below 0.001. In view of the utility of different cp regions, phylogenetic analyses were performed for the CCG data and the CDS sequences respectively.

GC Content Distribution and the Ka/Ks Ratios of Species Pairwise

The total and the third position GC content were compared between 41 species (belonging to Asparagaceae, Allioideae, Agapanthoideae, Asphodeloideae, Asphodelaceae, and Iridaceae). Lower GC contents were observed at the total nucleotides level (<38.5%) and the third codon positions (<36.0%) in most of Allium (Allioideae) species compared to species in non-Allioideae families (Figure 3 and Supplementary Table S6).

Open in a separate window

FIGURE 3

Changes in plastid GC content of all 41 taxa. This graph shows the total GC content (red bar and black line) and the third codon position GC content (blue bar and gray line) of each species.

The pairwise Ka/Ks ratios of each species pair were calculated (Figure 4), which provided information of selective pressure that acted on individual sequences. Much higher pairwise Ka/Ks ratios were observed in Allium (Allioideae) species pairs than non-Allioiseae species pairs (Figure 4 and Supplementary Table S7). In addition, high Ka/Ks ratios were also detected in other species (e.g., species in Hosta and Cordyline) (Figure 4 and Supplementary Table S7).

Open in a separate window

FIGURE 4

Pairwise Ka/Ks ratios in Allium (Allioideae) and other families. This heatmap shows pairwise Ka/Ks ratios between every sequence in the multigene nucleotide alignment. Allium (Allioideae) are shown on red branches. The scale factors associated with each value are shown on the right-hand side of the figure.

Repeat Sequences Variations, the Nucleotide Diversity, Codon Usage and Gene Loss

We detected numerous microsatellites (SSRs) in the six Allium cp genomes, ranging from 179 (A. maowenense) to 193 (A. chrysocephalum) (Supplementary Figure S1). The most abundant were mono-nucleotide repeats, where the number varied from 63 in A. chrysocephalum to 74 in A. rude, followed by tetra-nucleotides and di-nucleotide repeats, while the penta-nucleotide repeats were the least abundant (Supplementary Figure S1 and Supplementary Table S8). The overall length of the five categories of perfect SSRs ranged from 9 to 25 bp in the six Allium species (Supplementary Table S8).

The nucleotide diversity values in the LSC regions ranged from 0.0006 to 0.07823 with a mean value of 0.0310 (from 0.0035 to 0.0722 with the average value was 0.0465 in SSC regions), while the value was from 0.0000 to 0.0311 with a mean value of 0.0084 in the IRs regions (Supplementary Figure S2). Six genes with high nucleotide diversity were detected (>0.0700), these were trnK-UUU, matK, trnG-UCC, trnG-GCC, ndhF and rps15. Five genes (i.e., accD, clpP, rpl16, ccsA and ndhA) with nucleotide diversity more than 0.05500 were also detected.

The pattern of codon usage bias in the Allium (Allioideae) and non-Allioideae were investigated. We found that five parameters involved in codon usage bias were lower in Allium (Allioideae) species than non-Allioideae species, except the CAI that was lower in the family Asparagaceae (Figure 5).

Open in a separate window

FIGURE 5

The comparative analysis of codon usage bias in Allium (Allioideae) and other family species. (A) Labels representing each of families. (B) CAI (Codon adaptation index), (C) CBI (Codon bias index), (D) FOP (Frequency of optimal codons index), (E) GC3s (GC of synonymous codons in 3rd position), (F) ENC (Effective number of codons).

As shown in Figure 6, we found that the gene rps2 was lost in all Allium (Allioideae) species. In addition, four genes infA, rps16, psbZ, and ndhD were lost in four Allium species with different degree (A. sativum, A. macleanii, A. platyspathum, and A. victorialis). Gene cemA, infA, rps19, and ycf1 were lost in some species of Asparagaceae, and gene rpl32 and infA were missing in Aloe vera and A. maculate of Asphodeloideae.

Open in a separate window

FIGURE 6

Loss of chloroplast protein-coding genes in the phylogeny of all 41 taxa. Below is the phylogeny of all 41 species based on chloroplast genomes as shown in Figure 7. Different chloroplast regions were indicated at the left side. Allioideae are shown on red branches. IR, inverted repeat; LSC, large single-copy region; SSC, small single-copy region.

Characteristics of cp Genome and Phylogenetic Analysis

Complete chloroplast genome of the six species in section Daghestanica were newly sequenced in this study, and were deposited in GenBank (Supplementary Table S3). These plastid genomes are similar to previously published Allium plastomes in size, structure and gene content (Filyushin et al., 2016, 2018; Lee et al., 2017; Jin et al., 2018). The CCG data set had an aligned length of 192056 bp, within which 24727 parsimony-informative sites (PICs, 12.87 %) were detected. The SCG (single copy genes) possessed 54401 bp aligned nucleotides with 6755 PICs (12.41 %). CDS (coding DNA sequences) data set had an aligned length of 45954 bp nucleotides with 5464 PICs (11.89 %). Comparing these data sets with CCG, the percentages of PICs, SCG and CDS were reduced.

We reconstructed separate phylogenetic trees based on different methods: Bayesian and ML analyses on CCG. Bayesian and ML analyses recovered almost identical trees from each data set. There was strong support for the monophyly of each family were revealed based on CCG data (Figure 7). The topological structures from the SCG and CDS are similar to that from CCG, and all lineages possess high bootstrap values (Supplementary Figure S3). Agapanthus coddii from Agapanthoideae had strong support to be a sister to the Allioideae, and Asparagaceae was supported to be the sister of Agapanthoideae and Allioideae. For Chinese species in section Daghestanica, A. maowenense was closely clustered with A. herderianum, and A. chrysanthum is sister to the A. rude. All six species were closely clustered in one lineage (Figure 7 and Supplementary Figure S3).

Open in a separate window

FIGURE 7

Phylogenetic relationships of section Daghestanica species with related 35 species based on the whole cp genomes. (A) Tree constructed by Bayesian inference (BI) and maximum likelihood (ML) with the posterior probabilities of BI and the bootstrap values of ML above the branches, respectively. Represent maximum support in all two analyses. (B) Families of each species belong to, color of bar is consistent with the species’ color.

The Phylogenetic Analysis of Chinese Species in Section Daghestanica

The results of our phylogenetic analysis strongly support that Allium is monophyletic, which is in accordance with previous studies (Friesen et al., 2006; Nguyen et al., 2008; Li et al., 2010). We also found that Agapanthus coddii was closely related to Allium (Allioideae) (Figure 7 and Supplementary Figure S3). The relationships of Chinese species in section Daghestanica were resolved: the position of A. herderianum was confirmed, which showed a close relationship with A. maowenense, and differentiated early. A. chrysanthum was tightly clustered with A. rude, which is inconsistent with Li et al. (2010), who showed that A. chrysanthum was closely related to A. xichuanense, and A. rude was clustered with A. chrysocephalum. Our results may be more reliable, since we used the whole chloroplast genome (CCG) (Figure 7), SCG and CDS data (Supplementary Figure S3) respectively in phylogeny reconstruction compared to the rps16 fragment in study of Li et al. (2010), and our results were also supported by the morphological characteristics in Yu et al. (2018), which showed similar testa cells of A. chrysanthum and A. rude that are not parallel and irregular in long axis, and their filaments are longer than perianth segments (Figure 1). Although the A. herderianum has a close relationship with A. maowenense, the morphological characteristics of them are distinct different, and A. maowenense is most differentiated in morphology among the six species in terms of its perianth color, green midrib in perianth (Figure 1, ​,7).7). In addition, Yu et al. (2018) found that the testa cells of A. maowenense are parallel and irregular, which are obviously different from the other five species. The outer perianth segments of A. chrysocephalum and A. herderianum are boat-shaped, and their style are longer than perianths, these characteristics are distinct from the other four species (Figure 1). The leaves of A. chrysocephalum are flat falcate and in A. herderianum are semiterete and fistulose, easily differentiating the two species. However, we did not find a close relationship between A. chrysocephalum and A. herderianum (Figure 7). Therefore, our study uncovered a new relationship of Chinese species in Allium section Daghestanica using the whole cp genome data.

The Adaptation Evolution of Allium (Allioideae) Plastome

A previous study using 150 species cp genomes showed that the GC content of these species ranged from 19.5 to 42.1% (Smith, 2009). We found that Allioideae species typically exhibited lower GC content than other non-Allioideae families’ species (Figure 3). There are two reasons that may result in lower GC content of plastid DNA, firstly, a neutral mutation process such as AT-mutation pressure or AT-biased gene conversion, which will reduce the GC content (Howe et al., 2003; Kusumi and Tachida, 2005; Khakhlova and Bock, 2010), and secondly, that selection for translational efficiency may lead to the lack of G and C observed in plastid genomes (Morton, 1993, 1998; And and Voelker, 2003). We suggest that mutation pressure from the evolution process may be a crucial factor resulting in the low GC content of Allioideae. The Ka/Ks ratios of Allioideae species are exceptionally high compared with those observed within non-Allioideae families’ species and this may be an indication of an elevated mutation rate in Allioideae plastid DNA (Figure 4). Mutations during evolutionary processes leading to reduced GC content has also been found in mitochondrial DNA, nucleomorph DNA, and in the genomes of symbionts, parasites, and pathogenic bacteria (Ogata et al., 2001; And and Voelker, 2003; Lane et al., 2007; Smith and Lee, 2008).

Elevated pairwise Ka/Ks ratios were observed in Allioideae species pairs compared to non-Allioideae species pairs (Figure 4 and Supplementary Table S7). The elevated Ka/Ks ratios are unlikely to be explained by changes in codon preference since we did not obtain obvious codon usage bias in Allioideae species (Figure 5). It is possible that other factors (e.g., habitat environment and adaptation evolution) may have contributed to the elevated Ka/Ks ratios. Allioideae is a variable group that is spread widely across the Holarctic region from the dry subtropics to the boreal zone (Friesen et al., 2006; Li et al., 2010). Furthermore, species in Allioideae grow in various conditions from dry and well-drained soils to moist and organic soils, with most growing in sunny locations, and a number of species also grow in forests, or even in swamps or water (Block, 2010); The environment imposes stressful living conditions on Allioideae species and results in species divergence. The Ka/Ks ratios have been widely used to infer the evolutionary dynamics and identify adaptive signatures among species (Yang and Bielawski, 2000; Fay and Wu, 2003; Ai et al., 2015), and elevated Ka/Ks ratios indicate species may have undergone more selective forces (Hurst, 2002). Thus, the elevated Ka/Ks ratios observed throughout the Allioideae may suggest species in Allioideae undergo some selection pressure that is unknown.

Gene losses and gains are considered as important adaptive processes that greatly contribute to trait evolution (Hahn et al., 2007; Ding et al., 2012). In this study, we found that the gene rps2 was lost in all Allium (Allioideae) species, with infA, rps16, psbZ, ndhD, cemA, rps19, ycf1, and rpl32 were lost in some species of Allioideae, Asparagaceae and Asphodeloideae (Figure 6). Gene infA, which codes for translation initiation factor 1, was lost in an early ancestor of Fabales and Cucurbitales (Millen et al., 2001), and it was found as pseudogene in many genus (e.g., Albuca, Behnia, Camassia, and Echeandia) in study of McKain et al. (2016). The study of Steele et al. (2012) identified that gene rps16, rpl32, and rps19 were missing from various taxa throughout Asparagales, and these shared losses were suggested as the result of common ancestry. The gene rps2 was also identified as a pseudogene in Chlorophytum rhizopendulum (McKain et al., 2016). However, the mechanism underlying the loss of the rps2 gene in Allioideae was poorly understood. A previous study indicated that the product of the rps2 gene plays an important role in defense signal transduction (Bent et al., 1994). Although we detected the loss of the rps2 gene in Allioideae, we did not find a cause. Therefore, further studies are needed to examine whether specific factors were associated with the loss of the rps2 genes in the Allioideae.

We acknowledge Hao Li, Fu-Min Xie, and Xin Yang for their help in materials collection. We would like to thank Jun Wen, Juan Li, and Jiao Huang for their help in software use.