Comparative analysis and implications of the chloroplast genomes of three thistles (Carduus L., Asteraceae)

Taxon sampling, total DNA extraction, chloroplast genome assembly, and comparative analysis

Phylogenetic analysis of Carduus and related taxa

A total of 78 protein-coding regions were extracted from the complete cpDNA of the focal species and other related taxa (Table S1). Sequences were aligned using MUSCLE embedded in Geneious Prime (Edgar, 2004). We used jModeltest (Posada, 2008) to find the best model for the aligned DNA sequences; GTR + I + R was selected as the most suitable model and was used in Maximum Likelihood (ML) and Bayesian Inference (BI) analyses. The ML analysis was conducted with the IQ-tree web server (http://iqtree.cibiv.univie.ac.at), using 1,000 bootstrap replications to calculate branch support values (Trifinopoulos et al., 2016). We used MrBayes v3.2 (Ronquist et al., 2012) for BI analyses. The Markov chain Monte Carlo (MCMC) analysis was run for 1,000,000 generations, and a tree was assembled every 1000 generations. A 25% burn-in setting was used for summarising trees. Figtree v4.0 (http://tree.bio.ed.ac.uk/software/figtree/) was used to visualise phylogenetic trees. Other datasets, including whole chloroplast genomes (excluding one IR region), non-coding regions of cpDNA, and hotspot regions derived from the cpDNA of Carduus species, were used in phylogenetic analysis in addition to protein-coding regions (Table S1). Analytical procedures for these additional datasets were identical to those used for the protein coding regions; however, we used the TVM+I+G model for the whole chloroplast genome (excluding one inverted repeat [IR] region) and all non-coding regions, and the TVM+G model for the hotspot regions dataset.

SNP identification, primer design, and multiplex PCR

The complete matK gene, extracted from the cpDNA of the three focal species, was aligned using MUSCLE to identify SNPs (Edgar, 2004; Fig. S1). The selected SNP for C. crispus was then confirmed by aligning the available matK sequences of other Carduus species on NCBI to those of the focal species (Fig. S2). Based on SNP data, primer pairs were designed using Primer3 to distinguish C. crispus from other Carduus species (Untergasser et al., 2012). Primer sequences included matK_463F (5′-CATCTGGAAATCTTGGTTCAG-3′), matK_1162R (5′-GATGCCCCAATGCGTTACAA-3′), CD_SNP_F1 (5′-AATTCTTGCTTCAAAAGG GTCC- 3′), CD_SNP_R1 (5′-TTCCATTTATTCATCAA AAGATAC-3′), CD_SNP_F2 (5′-AATTCTTGCTTCAAAAGGGTCG-3′), and CD_SNP_R2 (5′-TTCCATTTATTCATCAAAAGATAG- 3′). The multiplex PCR of matK_463F, matK_1162R, CD_SNP_F1, and CD_SNP_R1 was designed to yield the 323 bp band for C. crispus, the 421 bp band for other Carduus, and the 700 bp band for all examined samples (Fig. S3). By contrast, the combination of matK_463F, matK_1162R, CD_SNP_F2, and CD_SNP_R2 yielded a 421 bp PCR product for C. crispus and a 323 bp band for other Carduus species (Fig. S3). Reactions were conducted in 25 µl solution consisting of 50 ng of template DNA, 2.5 µl of 10× reaction buffer, 0.5 U of E-taq DNA polymerase, 50 mM MgCl₂, and 5 mM dNTPs. Concentrations of outer primer pairs (matK_463F and matK_1162R) and inner primer pairs (CD_SNP_F1 and CD_SNP_R1, CD_SNP_F2 and CD_SNP_R2) were 0.75 pM and 0.5 pM, respectively. The PCR procedure consisted of 1 min at 94 °C, followed by denaturing for 1 min at 94 °C, annealing for 40 s at 55 °C, an extension stage of 50 s at 72 °C, and an additional extension of 7 min at 72 °C.

Comparative chloroplast genome analysis of the focal species

Differing numbers of reads were obtained from NGS data, resulting in varying cpDNA coverage rates among the three focal species (Table 2). Total cpDNA length differed among species, ranging from 152,342 bp to 152,617 bp, and included a large single copy (LSC), a small single copy (SSC), and two IR regions (Fig. 1). By contrast, all three species had identical numbers of protein coding (80), tRNA (30), and rRNA (4) genes (Table 2, Table S3). The IR-LSC and IR-SSC junctions were located in the rps19 and ycf1 coding regions, respectively, but were longer in C. acanthoides (rps19 = 61 bp and ycf1 = 568 bp) than in the other two species (rps19 = 60 bp and ycf1 = 565 bp). In addition, pairwise identity indicated that C. crispus is more similar to C. tenuiflorus (99.6%) than C. acanthoides (99.2%). Observations of nucleotide diversity indicated that 119 of 131 surveyed regions differ among the three focal species (Table S2, Fig. 2). Compared to coding regions, non-coding sequences had higher Pi values (Fig. 2). The highest Pi values were found in the psbC-trnS (0.0171) and psbH-petB (0.0161) regions. The highest value in coding regions was 0.00696 for ycf1 (Fig. 2). High nucleotide diversity regions (Pi values >0.008) included psbI-trnS_GCU, trnE_UUC-rpoB, trnR_UCU-trnG_UCC, psbC-trnS_UGA, trnT_UGU-trnL_UAA, psbT-psbN, petD-rpoA, and rpl16-rps3.

Features of cpDNA repeats

Analysis of SSRs yielded 43 SSRs in C. crispus, 40 in C. tenuiflorus, and 31 in C. acanthoides (Table S4). SSRs occupied the same position in all three species and were mostly located in non-coding regions. Although four types of SSR (i.e., mono-, di-, tri-, and tetra-nucleotides) were identified, most SSRs were mononucleotides composed of A and T nucleotides (Table S4). All 31 SSRs found in C. acanthoides were also present in the other two species (Table S4). By contrast, C. crispus had three unique SSRs and shared nine SSRs with C. tenuiflorus. There were no unique SSRs in C. acanthoides or specific shared SSRs between C. acanthoides and either C. crispus or C. tenuiflorus.

Among the focal species, 16 repeats were identified for both C. crispus and C. tenuiflorus, compared to 15 for C. acanthoides (Fig. 3A, Table S5). There were more repeats in coding regions than in non-coding areas, with the exception of C. acanthoides. Three types of repeats (i.e., forward, reverse, and palindrome) were identified in C. tenuiflorus and C. acanthoides; by contrast, only forward and palindrome repeats were found in C. crispus. Forward repeats were more abundant than reverse and palindrome repeats (Fig. 3B). Among recorded repeats, nine were common among all three species (Fig. 3C). Carduus acanthoides had four unique repeats, whereas C. tenuiflorus and C. crispus each had a single unique repeat. Carduus acanthoides shared one specific repeat with both C. crispus and C. tenuiflorus, whereas C. crispus and C. tenuiflorus shared five specific repeat regions (Fig. 3C).

Phylogenetic relationships between Carduus and related taxa

The ML and BI analyses, based on 78 protein-coding genes from Carduus and related taxa, yielded identical topologies (Fig. 4). In particular, both analyses confirmed the monophyly of Asteraceae subfamilies (i.e., Carduoideae, Chichorioideae, and Asteroideae). In contrast to the high support for Carduoideae and Cichorioideae (Bootstrap = 100/Posterior Probability = 1), low support was found for Asteroideae clades (Fig. 4). Notably, monophyly of the three Carduus species was not supported by either analysis. For example, C. acanthoides was sister to Silybum marianum, whereas C. crispus and C. tenuiflorus were sister to Cirsium arvense. Additional ML and BI analyses of full cpDNA sequences and non-coding regions suggested similar relationships (Fig. S4).

Multiplex PCR and specific markers for C. crispus

The results of multiplex PCR for the two groups of primer pairs yielded similar products, both of which were designed to identify C. crispus. In the first group, a 323 bp band was found in C. crispus, whereas a 421 bp band was identified in C. acanthoides and C. tenuiflorus (Fig. 5A). By contrast, the combination of matK_463F, matK_1162R, CD_SNP_F2, and CD_SNP_R2 yielded a longer PCR product for C. crispus in comparison with the other species (Fig. 5B). The designed primer pairs were specific to all C. crispus samples examined in this study (Figs. S5 and S6).

Conservatism of Carduus cpDNA

Chloroplast genomes are highly conserved in angiosperms with respect to gene content and order (Sugiura, 1992). This conservative tendency was observed in the three newly sequenced Carduus cpDNA genomes, compared to other Asteraceae members (Table 2). Other cpDNA sequences have revealed unique genomic events in Asteraceae. For example, the atpB gene, which encodes the CF1 ATPase beta subunit, is annotated as a pseudogene in Aster spathulifolius due to a deletion within the coding region (Choi and Park, 2015). Similarly, trnT_GGU was completely deleted or pseudogenised in the tribe Gnaphalieae (Lee et al., 2017). Duplication of trnF_GAA has been identified in Taraxacum (Salih et al., 2017). No comparable genomic events are present in Carduus or other members of subfamily Carduoideae (Table 2, Fig. S7). However, nucleotide diversity data pointed to potential regions for further study of phylogeny and population genetics, and the development of Carduus-specific molecular markers (Table S2, Fig. 2). The number of species we examined for this study was low relative to the approximately 32,000 known Asteraceae species. Therefore, additional studies that include the majority of Asteraceae species should be conducted to explore the overall evolutionary trends in the chloroplast genomes of this globally-distributed family.

Chloroplast genomes provide useful molecular data for reconstructing phylogeny, exploring biogeography, and estimating divergence time in angiosperm lineages (Do, Kim & Kim, 2014; Nguyen, Kim & Kim, 2015; Kim, Kim & Kim, 2016; Kim & Kim, 2018). Repetitive sequences in the chloroplast genome provide useful information for studying genomic rearrangement and phylogeny (Cavalier-Smith, 2002; Nie et al., 2012; Yi et al., 2013; Kim & Kim, 2018). In addition, existing repeats might result in the accumulation of new repeats in cpDNA (Asano et al., 2004). One of the crucial molecular data types in cpDNA is SSR sequences. Other studies have used SSRs to develop specific markers for different species and to study the genetic diversity of angiosperms (Ishikawa, Sakaguchi & Ito, 2016; Luo et al., 2016; Marochio et al., 2017; Han et al., 2018). In this study, although cpDNA sequences were highly conserved, the three Carduus species were found to have different numbers of SSRs (Table S4). While we did not develop SSR markers or conduct population studies of Carduus, the SSR information we obtained may be useful in future studies on Carduus species. In addition to the repeats shared among the three species, C. crispus had three unique repetitive sequences (Table S4), which may be useful in population studies, phylogenetic analyses, and the development of additional molecular markers.

Uncertain relationships among Carduus

Phylogenetic analyses of the Asteraceae have identified ambiguous relationships between Carduus, Cirsium, and Silybum (Fu et al., 2016; Panero, 2016; Arnelas et al., 2018); for example, ITS data suggests that Carduus is polyphyletic (Häffner & Hellwig, 1999). Although three coding regions (matK, rbcL, and ndhF) were used to reconstruct phylogenetic relationships, the position of Carduus remained unresolved (Fu et al., 2016). We used 78 protein-coding regions to clarify these relationships; however, the phylogeny of Carduus and related taxa remains unclear (Fig. 4). Specifically, C. acanthoides was found to be close to Silybum marianum whereas C. crispus and C. tenuiflorus form a clade with Cirsium arvense. While non-coding regions can be useful in reconstructing the phylogeny of lower taxa, we were unable to recover the monophyly of Carduus using data from non-coding regions, including the eight hotspot areas as well as the combined data from coding and non-coding regions (Fig. S4). These issues suggest the need for additional studies on the phylogeny of Carduus and other members of the subfamily Carduoideae using supplementary molecular data and morphology.

Implications of SNP data for developing molecular markers for Carduus

SNPs are useful in population studies due to its extremely abundant presence in the angiosperms genomes (Cui et al., 2017; Fischer et al., 2017; Pantoja et al., 2017), and are effective in phylogenetic analysis (Leaché & Oaks, 2017). In addition, various molecular markers have been developed for different angiosperm species based on SNP data from chloroplast genomes (Khlestkina & Salina, 2006; Wang et al., 2015a; Wang et al., 2015b; Hyun et al., 2019; Xia et al., 2019). We successfully developed a molecular marker, inferred from SNP data, to distinguish C. crispus from C. acanthoides and C. tenuiflorus (Fig. S3). Our marker demonstrates that nucleotide sequence variations can provide rapid molecular identification of C. crispus. We focused on C. crispus because it exhibits the characteristics of an invasive species (Dunn, 1976; Verloove, 2014; Jung et al., 2017), and may also have value for the treatment of obesity and cancer (Xie, Li & Jia, 2005; Davaakhuu, Sukhdolgor & Gereltu, 2010; Lee et al., 2011; Tunsag, Davaakhuu & Batsuren, 2011). Various DNA-based markers (i.e., inter-simple sequence repeats [ISSRs], sequence characterisation of amplified regions [SCARs], and SSRs) have been developed to authenticate medicinal plants to ensure safety and efficacy (Hao et al., 2010; Sarwat et al., 2012; Ganie et al., 2015; Ward, Gaskin & Wilson, 2008). Additionally, molecular data are useful for understanding invasion processes of alien plants (Ward et al., 2008). We developed a SNP-based molecular marker for C. crispus (Fig. 5, Figs. S5 and S6) that may be used to detect C. crispus invasions in their early stages, and to develop suitable management strategies. Our alignment results also identified specific SNPs for C. acanthoides and C. tenuiflorus, which may be used to create molecular markers for these species (Fig. S1). Although we only used the SNP in matK, use of the complete cpDNA sequences of Carduus will enable the mining of SNPs from other regions for developing molecular markers for C. crispus and related species.