A comparative analysis of the complete chloroplast genomes of three Chrysanthemum boreale strains

Ethics statement

The plant sample used in this study is neither endangered nor protected, and was collected from an area that was not privately owned or protected in any way. No specific permits were required to conduct this study.

Plant materials and sequencing

Two C. boreale strains with morphological differences were collected from different locations (Fig. S1) in the Republic of Korea and deposited at the National Agrobiodiversity Center, Rural Development Administration. The strain from Gongju-si, Chungcheongnam-do was labeled IT232531, and the one from Suwon-si, Gyeonggi-do was labeled IT301358. The total DNA was isolated from fresh leaves as previously described (Kim et al., 2006). The quality and quantity of DNA were examined using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and gel electrophoresis (in 0.8% agarose). Paired-end libraries of 350-bp insert size were constructed using TruSeq DNA PCR-Free kit (Illumina, San Diego, CA, USA) and sequenced with a 101-bp read length by Macrogen (Republic of Korea) using the HiSeq4000 (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. Another C. boreale strain, labeled 121002, was collected from Jeongeup-si, Jeollabuk-do (Hwang et al., 2013) and its cp genome was sequenced. Our group had previously submitted this cp genome sequence to NCBI with accession number MG913594 (Won, Jung & Kim, 2018).

Chloroplast genome assembly and annotation

The complete cp genome was assembled de novo (Kim et al., 2015a). Briefly, raw reads were trimmed using the Trimmomatic program (Bolger, Lohse & Usadel, 2014), assembled using the clc_assembler in the CLC Genomics Workbench v6.0 (CLC Bio, Denmark, Europe). Gaps were filled using Gap Closer (Luo et al., 2012). The resulting contigs were searched for cp-encoding contigs by BLASTN analysis against the cp genome of C. boreale strain 121002, and circularized. These were annotated using the online programs Dual Organellar GenoMe Annotator, cpGAVAS v.2.0 and BLAST (Shi et al., 2019; Wyman, Jansen & Boore, 2004). The structure of transfer RNA (tRNA) was predicted using the tRNAscan-SE 1.21 program using the default settings (Schattner, Brooks & Lowe, 2005). The circular genome map with structural features was generated using the OGDRAW v1.2 program (Lohse et al., 2013). The resulting cp genome sequences of strains IT232531 and IT301358 were deposited in NCBI under the IDs MN909052 and MN913565, respectively.

Chloroplast genome comparison

The cp genomes of the three C. boreale strains were compared using the mVISTA program in the Shuffle-LAGAN mode, using the annotation of strain 121002 as the reference (Frazer et al., 2004). The SNPs and indels in the cp genome were also recorded using DnaSP6.0 (Rozas et al., 2017) and manually verified from the sequence alignment by Clustal Omega (Sievers et al., 2011).

Characterization of repetitive sequences

Simple sequence repeats (SSRs) were discovered using the online web tool MISA (http://pgrc.ipk-gatersleben.de/misa/) with the following parameters: ten repetitions for mononucleotide motifs, eight for dinucleotide motifs, four for tri- and tetra-nucleotide motifs, and three for penta-and hexa-nucleotide motifs (Beier et al., 2017). Next, four different types of repeats, namely forward (F), palindromic (P), reverse (R) and complement (C) repeats were analyzed using the REPuter program (https://bibiserv.cebitec.uni-bielefeld.de/reputer) with a minimum repeat size of 30 bp and a Hamming distance of 3 (Kurtz et al., 2001). To reduce redundancy, IRb sequence was removed before analysis and repeats detected at the same position were merged into single repeat.

Phylogenetic analysis

The entire cp genomes and 77 protein-coding sequences shared in cp genomes of species belonging to tribe Anthemideae were used to reconstruct the phylogenetic relationships. Lactuca sativa L. was used as the outgroup. The species and the accession numbers of their cp genomes in NCBI are listed in Table S1. The nucleotide sequences were aligned using Clustal Omega (Sievers et al., 2011). Maximum likelihood (ML) analyses were conducted using the IQ-TREE web server (http://iqtree.cibiv.univie.ac.at) with the best-fit models determined by ModelFinder in the IQ-TREE package (Table S2) and 1,000 bootstrap replicates (Hoang et al., 2018; Kalyaanamoorthy et al., 2017; Nguyen et al., 2015). Bayesian inferences (BI) were performed with MrBayes v. 3.2.7 (Ronquist et al., 2012) and the nucleotide substitution models determined by ModelTest-NG (Darriba et al., 2019) (Table S2). The Markov chain Monte Carlo algorithms were run for 10 million generations and sampled every 1,000 generations. The first 25% of trees were discarded as burn-in and the remaining trees were used to build a majority-rule consensus tree with posterior probability values for each node. The stationary was considered to be reached when the average standard deviation of split frequencies remained below 0.01. The phylogenetic tree was visualized with FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Characterization of chloroplast genomes

We used NGS techniques to generate approximately 30.2 Gb and 34.9 Gb of raw reads from strains IT232531 and IT301358, respectively. We assembled de novo the complete cp genomes of sizes 151,098 bp for IT232531 and 151,010 bp for IT301358 (Table 1). For comparison, we included the previously reported cp genome of C. boreale strain 121002, which was 151,012 bp in size (Won, Jung & Kim, 2018). All the three C. boreale strains had a typical quadripartite structure of cp genomes with an LSC, an SSC, and a pair of IR regions (Fig. 1). The length of the LSC region was 82,817 bp, 82,880 bp and 82,788 bp for the strains 121002, IT232531 and IT301358, respectively. The SSC region measured 18,281 bp, 18,312 bp and 18,310 bp in the strains 121002, IT232531 and IT301358, respectively. The strains were comparable in terms of the length of the IR regions and the GC content of the LSC, SSC, IR regions and the complete genome (Table 1). The IR regions had a higher GC content than the LSC and SSC regions due to the presence of GC-rich ribosomal RNA (rRNA) genes and tRNA genes in these regions.

The cp genomes of all the strains comprised 113 unique genes. These included 80 protein-coding genes, 29 tRNA genes and four rRNA genes (Table 2). Each strain contained 61 protein-coding genes and 21 tRNA genes in the LSC region and 11 protein-coding genes and one tRNA gene in the SSC region (Fig. 1). Three genes (rps12, rps19 and ycf1) were distributed in both single copy and IR regions. The IR regions contained seven protein-coding genes, seven tRNA genes and four rRNA genes each. Because the IR regions are duplicates of each other, all genes in these regions were also duplicated.

The cp genomes of C. boreale included 16 intron-containing genes (Table 3). The genes ycf3 and clpP had two introns each, while all other genes contained a single intron. Nine of the introns were identical in length, whereas seven other introns differed in length between 1 bp and 24 bp. The intron of the trnK-UUU gene was largest (2,560–2,575 bp) in all the strains and its pairwise length differed between the strains by 7–15 bp. The intron of the ndhA gene in IT232531 was 24 bp and 6 bp longer than that in strains 121002 and IT301358, respectively. In each strain, the rps12 gene was trans-spliced, with the 5′ end exon located in the LSC region and the duplicated 3′ end exon located in both the IR regions, as previously reported in other plants (Thode & Lohmann, 2019).

Given that the cp genome of C. boreale strain 121002 was obtained using PacBio’s long reads (Won, Jung & Kim, 2018), we repeated the cp genome assembly of 121002 using Illumina’s short reads as conducted for other C. boreale strains. The sequence comparison between two cp genomes obtained using long reads and short reads revealed that there was no SNP detected. Instead, indels were observed at four genomic regions and all of them were associated with homopolymers. Three indels were located in intergenic spacers (IGSs), trnE-UUC_rpoB (18 thymines in the long-read assemble vs. 17 thymines in the short-read assemble) and psaA_ycf3 (16 vs. 15 adenines), and the intron of rpl16 (8 vs. 9 cytosines), which didn’t change the protein sequences. However, the coding region of ycf1 possessed one indel (13 adenines vs. 14 adenines) (Fig. S2), which resulted in 1,036 amino acids (aa) in the original data due to the premature stop codon. However, one-bp insertion generated the ycf1 protein of 1,668 aa, which was more consistent with the other C. boreale strains (1,672 aa in IT232531 and 1,673 aa in IT301358). While we used the original cp sequence deposited in NCBI for analyses, in case of ycf1, we used the newly obtained sequences.

Variation in chloroplast genomes

The mVISTA-based identity plot indicated conservation in DNA sequence and gene synteny across the whole cp genome, and revealed the regions with increased genetic variation (Fig. 2). The gene number, order and orientation were conserved. There was higher genetic variability in the single copy (LSC and SSC) regions than in the IR regions, and in non-coding regions than in coding regions. Highly diverged regions included the IGSs, trnK-UUU_rps16, trnS-GCU_trnC-GCA, trnR-UCU_trnT-GGU, rps4_trnL-UAA, ndhC_trnV-UAC, psbE_petL, rps16_rps3, and trnL-UAG_rpl32 and the introns of trnK-UUU, rps16 and ndhA (Fig. 2). We detected a total of 298 SNPs (Table S3). The LSC region contained a majority of the SNPs (204, accounting for 68.5% of the SNPs), followed by the SSC region (75, 25.2%), and the IR regions (19, 6.4%). The SNPs were more abundant in non-coding regions: 141 were located in intergenic regions, 46 in introns, and 111 in coding regions. The ycf1 gene contained the largest number of substitutions (25 SNPs), followed by the trnK-UUU intron (18 SNPs), rpoC2 (12 SNPs) and the ycf1_rps15 IGS (11 SNPs).

We detected a total of 106 indels (Table S4): 81 in the LSC, 19 in the SSC, and six in the IR regions. A total of 86 and 17 indels were located in IGS and introns, respectively, whereas three were contained in coding regions. The ndhC_trnV-UAC spacer had five indels, while the introns of trnK-UUU and ndhA, and the spacers of psaA_ycf3 and psbE_petL contained four indels each. The psbE_petL IGS included the two largest indels (54 bp and 36 bp) in the cp genome. The trnK-UUU intron was the longest in the genome, and one of the most variable regions, comprising both SNPs and indels (Fig. S3). The 5-bp deletion at the end of the protein-coding gene rpoC2 in strain IT232531 generated a protein that was longer by two amino-acids. In the ycf1 gene, the 3-bp insertion in strain IT301358 did not change the protein’s translational frame.

We investigated the position of genes at the junction regions (LSC/IRa, IRa/SSC, SSC/IRb and IRb/LSC; Fig. 3). At the LSC/IRa junction, C. boreale possessed rps19 with 220 bp in LSC and 59 bp in IRa. The IRa/SSC junction contained the functional ycf1, while the SSC/IRb possessed the duplicated partial copy, pseudogene ycf1 (Ψycf1) and ndhF. At the IRb/LSC junction, rpl2 and trnH-GUG were located within the distance of 122–124 bp from each other.

Repeat analysis

We investigated the distribution of SSRs that were 1–6 bp in length in the C. boreale cp genomes. We recorded a total of 47, 43 and 50 SSR motifs in 121002, IT232531 and IT301358, respectively (Table S5). Mononucleotide repetition was most prevalent in each cp genome, followed by tri-, penta-and tetra-nucleotide repetition (Fig. 4A). We did not detect di-or hexa-nucleotide SSRs. In terms of sequence context, there were more adenine and thymine residues than cytosine and guanine residues (Fig. 4A). Intergenic and intronic regions contained more SSRs than coding regions, with 41, 37 and 42 instances of SSR occurrence in the non-coding regions in the 121002, IT232531 and IT301358 strains, respectively (Fig. 4B). Most of the SSRs were located in the LSC region followed by those in the IR region (Fig. 4C).

We detected four different types of long dispersed repeats (LDRs), namely forward (F), palindromic (P), reverse (R) and complement (C) repeats, each with a motif length longer than 30 bp. We identified a total of 19 (12F, 5P, 2R), 18 (11F, 7P) and 16 (9F, 5P, 1R, 1C) repeats in the cp genomes of strains 121002, IT232531 and IT301358, respectively (Fig. 5A; Table S6). F and P repeats were more abundant than C and R repeats. Repeat units of 30–34 bp were the most common, whereas repeat units longer than 40 bp occurred less frequently (Fig. 5A). More LDRs were located in non-coding regions (IGS and introns) than in coding regions (Fig. 5B). Among the protein-coding genes, LDRs were detected in the psaA, psaB and ycf2 in all three C. boreale strains (Table S6). Most LDRs were present in LSC region compared to IR and SSC regions, while some LDRs were shared among LSC, IR and SSC regions (Fig. 5C).

Phylogenetic analysis

The phylogenetic trees were constructed based on complete cp genome sequences and 77 protein-coding genes that were common to the three C. boreale strains, the 17 other species of the tribe Anthemideae (Asteroideae, Astereaceae), and the outgroup species, L. sativa (Cichorieae, Cichorioideae, Asteraceae). The multiple alignment of complete cp genomes contained 158,397 nucleotide sites in which 11,844 were variable and 3,371 were parsimony informative. The multiple alignment of protein-coding sequences possessed 62,965 nucleotide sites in which 3,065 were variable and 900 were parsimony informative. In each cp sequences, both ML and BI trees revealed similar topologies but minor difference within Chrysanthemum species (Fig. 6; Fig. S4). Two datasets also resulted in similar phylogenetic relationship. All Chrysanthemum sequences were grouped into a single clade together with Opisthopappus taihangensis (Ling) C.Shih with high bootstrap support and Bayesian inference (Fig. 6). Three C. boreale strains were all non-monophyletic, which was also observed in two C. morifolium analyzed. Additionally, Artemisia species were clustered into two clades. Among them, seven species formed a monophyletic group, and other three were located in another clade and were closer to the Chrysanthemum clade.