Comparison of intraspecific, interspecific and intergeneric chloroplast diversity in Cycads.

Cycads are among the most threatened plant species. Increasing the availability of genomic information by adding whole chloroplast data is a fundamental step in supporting phylogenetic studies and conservation efforts. Here, we assemble a dataset encompassing three taxonomic levels in cycads, including ten genera, three species in the genus Cycas and two individuals of C. debaoensis. Repeated sequences, SSRs and variations of the chloroplast were analyzed at the intraspecific, interspecific and intergeneric scale, and using our sequence data, we reconstruct a phylogenomic tree for cycads. The chloroplast was 162,094 bp in length, with 133 genes annotated, including 87 protein-coding, 37 tRNA and 8 rRNA genes. We found 7 repeated sequences and 39 SSRs. Seven loci showed promising levels of variations for application in DNA-barcoding. The chloroplast phylogeny confirmed the division of Cycadales in two suborders, each of them being monophyletic, revealing a contradiction with the current family circumscription and its evolution. Finally, 10 intraspecific SNPs were found. Our results showed that despite the extremely restricted distribution range of C. debaoensis, using complete chloroplast data is useful not only in intraspecific studies, but also to improve our understanding of cycad evolution and in defining conservation strategies for this emblematic group.

Cycads are iconic relict species, regarded as “living fossils” because of their recognizable intermediate morphological traits between angiosperms and gymnosperms1. Cycads dominated the Mesozoic but their origin can be dated to the late Paleozoic (~265–290 Ma)23. However, molecular dating studies indicate that living cycad species could be not much older than ~12 Ma, rejecting both the hypothesized role of dinosaurs in generating extant diversity and the use of “living fossils” to describe current cycad species4. Cycads are distributed in tropical and subtropical regions of Africa, Asia, Oceania and America5. Ten genera and 344 species are currently accepted6, with Cycas containing roughly 40% of the species in the Near Threatened and Vulnerable categories in the IUCN red list7.

The genus Cycas L., in the monotypic family Cycadaceae, is the oldest genus of cycads, holding about 113 species5. More than 20 species are found in China58, with most of them endemic. Cycas debaoensis Y. C. Zhong & C. J. Chen, a critically endangered cycad species endemic to southwest China9, only occurs in 11 small populations near the border of Guangxi province and Yunnan province10. Previous studies have assessed genetic diversity in C. debaoensis using inter simple sequence repeat (ISSR) markers or nuclear microsatellites, showing limited gene flow among populations and low within-population diversity101112.

Chloroplasts (cps) are present in photosynthetically active green tissues and generally develop from proplastids in meristems or etioplasts after illumination of dark-grown tissues, and display a conserved structure of two inverted repeats (IR) separated by small (SSC) and large (LSC) single-copy regions13. Due to their natural abundance in plant cells (≈3–5% of the cell DNA content comparing to nuclear DNA)14, cp sequences are a versatile tool for plant identification (DNA-barcoding) and evolutionary studies. They have been used at small and large temporal scales in plants1516. The use of cps is a very powerful tool to reconstruct plant phylogenies and infer historical biogeographic patterns of diversification1718. However, only a limited number of regions in the chloroplast genome have been used to address evolutionary, taxonomic and biodiversity questions in Cycas19. With the rapid development of Next Generation Sequencing (NGS), it is now feasible to obtain the entire sequence of the chloroplast using a genome skimming approach, resulting in high resolution phylogenies and allowing for estimations of timing of historical diversification, biodiversity and extent of genomic divergence142021.

It is well known that genetic diversity can greatly vary between taxa, due to either different intrinsic characteristics (e.g. reproductive system, genome size and organization) or to extrinsic features (e.g. endemic vs. widespread species, young vs. old species)171822. In addition, the hypothesis of a linear accumulation of mutations in sequences across time (i.e. a molecular clock) has been refuted in many groups232425.

In this study, we analyse the C. debaoensis chloroplast as a reference together with molecular data of other Cycas species to identify potential DNA-barcode loci, and compared generic-level chloroplast features in cycads, to highlight the evolutionary history of this group. We also compared the chloroplast features of two individuals of C. debaoensis to provide new resources for marker development in this endangered species.

Results and Discussion

Genome size and features

Using genome skimming and reference-guided assembly, we reconstructed the 162,094 bp long chloroplast genome C. debaoensis_Jiang_DB-2015. The complete cp genome was submitted to GenBank under accession number KU743927. It was 2 bp longer than C. debaoensis (KM459003) due to two 1 bp indels in the LSC (Fig. 1). The two C. debaoensis individuals exhibited the typical composition of LSC, SSC regions and two IR copies of 88,854 bp, 23,088 bp, and 25,076 bp (Table 1). The overall GC content of C. debaoensis_Jiang_DB-2015 was 39.4%, and 38.7%, 36.6% and 42.0% in the LSC, SSC and IR regions, respectively. These values are similar to C. debaoensis KM459003. C. debaoensis (KU743927), C. debaoensis (KM459003) and C. revoluta showed similar GC content (39.4%), slightly lower than C. taitungensis (39.5%) (Table 1). In total, 133 genes were annotated, including 87 protein coding genes, 37 tRNA genes and 8 rRNA genes in C. debaoensis_Jiang_DB-2015, while 156 and 169 genes were annotated in C. revoluta and C. taitungensis, respectively. Twelve genes (atpF, rpoC1, rpl2, clpP, ycf3, trnG-UCC, trnK-UUU, trnL-UAA in LSC; ndhA locates in SSC; ndhB, trnA-UGC, trnI-GAU in the IRs regions) contain 1–2 introns, respectively; while fourteen genes (3 protein coding 7 tRNA and 4 rRNA) were duplicated in the IR regions (Table 2). The tufA gene, found in gymnosperms and hornworts and inherited from green algae is coding for a nonfunctional protein synthesis elongation factor (723-bp long in C. taitungensis and Ginkgo)17. Interestingly, this gene was 1 bp longer in C. debaoensis_Jiang_DB-2015 than in other Cycas (Table 2).

Figure 1

Circular gene map of the plastid genome of Cycas debaoensis.

Genes drawn within the circle are transcribed clockwise, while those drawn outside are transcribed counter clockwise. Genes are colour-coded according to their functional groups. Inner circle: GC content.

Table 1

Characteristics of the complete chloroplasts used in the study.

Species	Length	LSC	SSC	IR	GC%				Coding	tRNA	rRNA
					Total	LSC	SSC	IR
C. debaoensis KU743927	162,094	88,854	23,088	25,076	39.4	38.7	36.6	42.0	87	37	8
C. debaoensis	162,092	88,852	23,088	25,076	39.4	38.7	36.5	42.0	87	37	8
C. revoluta	162,489	88,977	23,376	25,068	39.4	38.6	38.9	42.0	109	39	8
C. taitungensis	163,403	90,216	23,039	25,074	39.5	38.8	36.5	42.0	122	38	8

The Length, Large Single Copy (LSC) region; Small Single Copy (SSC) region; Inverted Repeats (IR), the GC content of each region, and the number of coding, tRNA and rRNA loci are shown. Data from GenBank for C. debaoensis (KM459003), C. taitungensis17 (NC_009618) and C. revoluta (NC_020319).

Table 2

Genes present in the C.debaoensis (KU743927) chloroplast genome.

	Group of genes	Gene names
1	Photosystem I	psaA, psaB, psaC, psaI, psaJ
2	Photosystem II	psbA, psbB, psbC, psbD, psbE, psbF, psbH, psbI, psbJ, psbK, psbL, psbM, psbN, psbT, psbZ
3	Cytochrome b/f complex	petA, petB, petD, petG, petL, petN
4	ATP synthase	atpA, atpB, atpE, atpFa, atpH, atpI
5	NADH dehydrogenase	ndhAa, ndhBa (×2), ndhC, ndhD, ndhE, ndhF, ndhG, ndhH, ndhI, ndhJ, ndhK
6	RubisCO large subunit	rbcL
7	RNA polymerase	rpoA, rpoB, rpoC1a, rpoC2
8	Ribosomal proteins (SSC)	rps2, rps3, rps4, rps7 (×2), rps8, rps11, rps12 (×2), rps14, rps15, rps16, rps18, rps19
9	Ribosomal proteins (LSC)	rpl2 a, rpl14, rpl16, rpl20, rpl22, rpl23, rpl32, rpl33, rpl36
10	Other genes	accD, ccsA, cemA, chlB, chlL, chlN, clpPb, infA, matK, tufa
11	Proteins of unknown function	ycf1, ycf2 (×2), ycf3b, ycf4, ycf12
12	Ribosomal RNAs	rrn4.5 (×2), rrn5 (×2), rrn16 (×2), rrn23 (×2)
13	Transfer RNAs	trnA-UGCa (×2), trnC-GCA, trnD-GUC, trnE-UUC, trnF-GAA, trnfM-CAU, trnG-GCC, trnG-UCCa, trnH-GUG (×2), trnI-CAU, trnI-GAUa (×2), trnK-UUUa, trnL-CAA (×2), trnL-UAAa, trnL-UAG, trnM-CAU, trnN-GUU (×2), trnP-GGG, trnP-UGG, trnQ-UUG, trnR-ACG (×2), trnR-CCG, trnR-UCU, trnS-GCU, trnS-GGA, trnS-UGA, trnT-UGU, trnV-GAC (×2), trnW-CCA, trnY-GUA

(×2) Two gene copies in the IRs.

aGene containing one intron.

bGene containing two introns.

Repeat and SSR analysis

Using REPuter, seven repeats were found in the chloroplast of C. debaoensis_Jiang_DB-2015, which were three forward (F) and four palindrome (P), with no reverse and complement repeats discovered (Table 3). The repeats were mainly distributed in the intergenic spacers of transfer RNA genes, some of them being located in transfer RNA itself (Table 3). Interspecific comparison and analysis in broader Cycadaceae showed that C. revoluta had the highest number of repeats (24), while C. debaoensis contained the fewest (7), and C. taitungensis contained an intermediate number of repeats (16) (Supplementary Table S1, Fig. 2). In contrast, the comparison of simple sequence repeats (SSRs) revealed a relative conservatism in their numbers, with congeneric species showing similarities in both numbers and spatial patterns of SSRs occurrence (Fig. 3). This is of particular value in Cycadaceae, where species are usually scarcely distributed, and in which diagnostic morphological characters are often poorly defined or absent. These generic molecular biomarkers have the potential to provide useful diagnostic data in redefining complex paraphyletic and polyphyletic species groups in the family182226, aiding directly in conservation efforts.

Table 3

Repeat sequences and their distribution found by REPuter in the C. debaoensis (KU743927) chloroplast genome.

	Repeat1 start (location)	Repeat2 start (location)	Size (bp)	Type	Region
1	88,852 (trnI-CAU)	137,037 (IGS chlL-trnN-GUU)	20,257	P	LSC; IRb
2	28,396 (IGS ropB-trnC-GCA)	28,470 (IGS ropB-trnC-GCA)	39	F	LSC
3	55,293 (trnM-CAU)	55,340 (IGS trnM-CAU-atpE)	32	P	LSC
4	48,603 (IGS trnF-GAA-ndhJ)	48,629 (IGS trnF-GAA-ndhJ)	31	F	LSC
5	113,680 (IGS trnN-GUU-ndhF)	113,680 (IGS trnN-GUU-ndhF)	30	P	IRa
6	113,680 (IGS trnN-GUU-ndhF)	137,236 (IGS chlL-trnN-GUU)	30	F	IRa; IRb
7	137,236 (IGS chlL-trnN-GUU)	137,236 (IGS chlL-trnN-GUU)	30	P	IRb

IGS: Intergenic spacer.

Figure 2

Repeat sequences in four chloroplast genomes of Cycas.

REPuter was used to identify repeat sequences with length ≥30 bp and sequence identity ≥90% in the chloroplast genomes. F and P indicate the repeat type F (forward) and P (palindrome), respectively. Repeats with different lengths are indicated in different patterns.

Figure 3

Number of simple sequence repeats in four chloroplast genomes of Cycas, classified by repeat type.

mono-: mononucleotide SSRs; di-: dinucleotide SSRs.

In other groups, the presence and nature of repeats have been shown to be of great value in evolutionary and population analyses2728. Microsatellites (SSRs) are useful markers for population genetics, conservation of endangered species and species delineation22293031, as previously highlighted for C. debaoensis1012. There were 39 SSRs in the chloroplast of C. debaoensis, 34 (87%) and 5 (13%) mono- and di-nucleotides SSRs, respectively (Table 4). These SSRs were mainly distributed in the IGS region (29; 74%), and the other 26% were distributed in CDS genes (Table 4). C. revoluta and C. taitungensis had similar SSRs numbers and locations than C. debaoensis, most of them being mononucleotides and distributed in the IGS (Fig. 3). Interestingly, C. revoluta lacked some SSR patterns (G and GA mono- and di- nucleotides, respectively) (Table S2). These diagnostic SSRs can be used in combination with nuclear SSRs developed in the genus for Cycadaceae conservation or reintroduction, species biodiversity assessments and phylogenetic studies in native or introduced areas1230313233.

Table 4

Simple sequence repeats in the C. debaoensis (KU743927) chloroplast genome.

No.	SSR type	SSR	size	start	SSR-containing region
2	p1	A	20	8307	IGS (trnQ-UUG-psbK); LSC
	p1	A	20	73073	clpP; LSC
1	p1	A	19	123912	IGS (ndhG-ndhI); SSC
1	p1	A	14	32963	IGS (trnE-UUC-trnH-GUG); LSC
1	p1	A	13	68694	IGS (tufA-trnH-GUG); LSC
2	p1	A	11	101911	IGS (rps7-rps12); IRa
	p1	A	11	148078	IGS (trnV-GAC-rps12); IRb
2	p1	A	10	1872	trnK-UUU; LSC
	p1	A	10	11723	IGS (trnG-UCC-trnG-UCC); LSC
1	p1	T	19	15487	IGS (atpH-atpI); LSC
1	p1	T	18	10133	IGS (trnS-GCU-ycf12); LSC
2	p1	T	16	57648	IGS (atpB-rbcL); LSC
	p1	T	16	84489	rpl16; LSC
3	p1	T	15	54643	IGS (trnM-CAU-ndhC); LSC
	p1	T	15	84795	rpl16; LSC
	p1	T	15	134345	IGS (ycf1-chlN); SSC
1	p1	T	13	74064	clpP; LSC
2	p1	T	12	5951	rps16; LSC
	p1	T	12	118369	IGS (rpl32-trnP-GGG); SSC
4	p1	T	11	83227	rps8; LSC
	p1	T	11	88618	IGS (rpl23-trnI-CAU); LSC
	p1	T	11	102859	IGS (rps12-trnV-GAC); IRa
	p1	T	11	149026	IGS (rps7-rps12);IRb
3	p1	T	10	63350	IGS (ycf4-cemA); LSC
	p1	T	10	69919	IGS (trnP-UGG-psaJ); LSC
	p1	T	10	125905	ndhA; SSC
1	p1	G	14	70451	IGS (psaJ-trnH-GUG); LSC
1	p1	G	11	52764	IGS (ndhC-trnH-GUG); LSC
2	p1	G	10	70921	IGS (rpl33-trnH-GUG); LSC
	p1	G	10	149995	IGS (rps7-ndhB); IRb
1	p1	C	14	44747	ycf3; LSC
1	p1	C	11	5229	rps16; LSC
2	p1	C	10	17681	IGS (rps2-rpoC2); LSC
	p1	C	10	100943	IGS (ndhB-rps7); IRa
3	p2	(TA)14	28	29832	IGS (petN-psbM); LSC
	p2	(TA)9	18	1453	IGS (psbA-trnK-UUU); LSC
	p2	(TA)6	12	15670	IGS (atpH-atpI); LSC
1	p2	(GA)6	12	68377	IGS (psbE-petL); LSC
1	p2	(AT)6	12	15808	IGS (atpH-atpI); LSC

Cycads phylogenetic reconstruction and comparison

In the maximum likelihood (ML) phylogenetic tree, all but two nodes were highly supported (bootstrap support ≥95), with the accessions of C. debaoensis closely related to the other Cycas species (Fig. 4A). Cycas spp. diverged first in the Cycadales, followed by Dioon, a clade containing Zamia, Ceratozamia and Stangeria. Bowenia diverged from the remaining Zamiaceae with relatively high support (bootstrap support 80%), Macrozamia being as sister to a clade containing both Encephalartos and Lepidozamia (Fig. 4A). This chloroplast phylogeny confirms the division of Cycadales into two suborders, each of them being monophyletic in our analyses (Fig. 4B), but contradict the current family delimitations34, with the family Stangeriaceae being polyphyletic with high support, in agreement with the most recent phylogenetic work in cycads35. However, we found each genus in Stangeriaceae grouping with one of the clade in Zamiaceae, contrary to Salas-Leiva et al.35, in which Bowenia diverged in a basal position for all cycads but Cycas and Dioon, and Stangeria group with Zamia (Microcycas being absent from our dataset). C. debaoensis diverged in a basal position of the genus, with C. revoluta and C. taitungensis being closely related. In addition, the branch leading to Cycas is longer than the branches leading to the other genera, in accordance with the hypothesis of a recent diversification during the Miocene19. These results demonstrate the suitability and efficiency of using complete chloroplast sequences in reconstructing the evolutionary history of cycads, as previously demonstrated for other groups3637.

Figure 4

Maximum Likelihood phylogenetic tree of the available chloroplast sequences in GenBank for Cycadales, plus the chloroplast sequences of Cupressus sempervirens, Taxus mairei, Pinus strobus, Araucaria heterophylla and Ginkgo biloba as outgroups (A). For readability and better understanding of the branches lengths, a zoom on the Cycadales family is shown (B).

Although the clear delineation of the genera in cycads are mostly due to the lengths of the sequences provided by complete chloroplast sequencing, the variability was unevenly distributed along the circular molecule (Fig. 5). Indeed, the ribosomal RNA genes as well as the region between psbA and rpoC1 genes (0–25 kb) and between ndhC and rpl20 (50–72 kb) of the cp sequences exhibited relatively few variations among cycads (Fig. 5). A cluster of four ndh genes (120–124.5 kb) appeared to be strikingly conserved. Overall, the level of variation increased with taxonomic distance, meaning that regions showing polymorphisms among Cycas species, exhibited higher variation among different cycad genera. However, some notable well-defined (<500 bp long) polymorphic regions departed from this assertion in ycf1 and ycf2 IRa. Indeed, the polymorphisms in these regions were higher among Cycas than among cycads as a whole. Interestingly, three SNPs were located in other regions (trnL-trnF, clpP intron 2 and ycf2 IRb) that exhibited a continuous increase in polymorphism levels across the family. ycf1 has been recently proposed as a barcode locus38, despite it not being present in all genera39, and was identified with clpP among the most variable loci in Parthenium spp.40. Finally, trnL-trnF was previously used in Cycas phylogenetic studies41, but also in species identification of trees and ferns4243. Here, we stress the need to further assess these loci as potential more informative substitutes to the official barcode loci44.

Figure 5

mVISTA percent identity plot of available cycad chloroplasts, using C. debaoensis as a reference.

Vertical scale indicates the percentage of identity ranging from 50% to 100%. Coding regions are in blue and noncoding regions are in pink. Cladogram redrawn from Fig. 5B; branch lengths are not representative of evolutionary changes; bootstrap support is indicated on the nodes. Black arrows indicate intraspecific variations in C. debaoensis.

Comparative interspecific chloroplast genomic analysis

Focusing on the three Cycas species available in GenBank, the mVISTA results showed that the four chloroplasts were highly conserved; however, the coding regions appeared to be globally more variable than the non-coding regions (Fig. 5 and Suppl. Fig. 1). Furthermore, the coding regions, e.g. rpoB, psbC, clpP (intron), ycf1, and ycf2; psbA-trnH and trnL-trnF intergenic spacers showed promising levels of variations for further development in applications such as DNA-barcoding or phylogenetic reconstruction.

Cycas have been considered as a difficult group for DNA-barcoding45. Previously, it was reported that the psbA-trnH spacer was highly variable in Cycadales except Cycas45, but trnL-trnF was used in Cycas for phylogenetic studies41. Although rbcL + matK were chosen as a two-locus DNA-barcode for their universality and efficacy in land plant44, it was not variable enough in Cycas (Fig. 5 and Suppl. Fig. 1). rpoB, rpoC1, and non-coding (e.g. atpF-atpH, psbK-psbI, and rpl32-trnL) regions have been shown to be variable enough at higher taxonomic levels46, but also within cycads4748. In light of this, psbC, clpP (intron), ycf1, and ycf2 should be considered as candidates for future phylogenetic studies in Cycas.

Intraspecific comparison

Comparing the two individuals of C. debaoensis, we found 10 SNPs and 1 indel, plus one “N” position in our data that didn’t allow us to confirm its status (Table 5). Their genomic positions are indicated in Figs 1 and 5. Six and four SNPs were located in IGS and coding regions, respectively, and one indel in the clpC intron. The genetic distance between the two individuals was 0.005553% (i.e. ≈1 SNPs/indels per 20 kb). This result is consistent with previous studies showing low within-population diversity in addition to limited gene flow among populations in C. debaoensis1011. Whittall et al.49 found a comparable level of intraspecific divergence in pines, irrespective of the rarity of the considered species. However, in the pest species Jacobaea vulgaris (Asteraceae), the intraspecific divergence was four times higher50, perhaps due to it’s short generation times as opposed to those prevailing in slow growing and long-lived trees. Further studies are still needed to determine whether or not intraspecific genetic diversity is linked to geographic ranges or the intrinsic characteristics of the taxonomic group.

Table 5

Intraspecific SNPs between two individuals of C. debaoensis (KU743927).

Position	Type	Nucleotide	location	location type
8432		A/T	trnQ-psbK	IGS
9889		C/T	trnS-ycf12	IGS
33707		G/T	psbD	coding
44878		C/G	ycf3	exon 1
48140		A/C	trnL-trnF	IGS
60132		G/T	trnR-accD	IGS
62218		A/C	psaI-ycf4	IGS
73073		A/.	clpP	intron 2
86731		C/T	rps19-rpl2	IGS
90307		G/T	ycf2	coding (IRa)
160640		A/C	ycf2	coding (IRb)

Conclusions

Comparing genomic diversity at different taxonomical, but also spatial and temporal scales, we were able to reconstruct a robust phylogeny for cycads, and to identify regions showing promising levels of variation at three levels (familial, generic and intraspecific rank). These regions can provide useful and alternative loci for species identification and population-based studies for conservation, ecology and evolution. Despite their restricted geographic ranges, we showed that several, potentially diagnostic intraspecific variations can be found in the chloroplasts of different individuals C. debaoensis, including 10 SNPs and 1 indel in as of yet unstudied regions. Comparing results from the three scales, four regions appeared to be variable at the three considered taxonomic scales, namely ycf3, clpP, psbD and the trnL-trnF IGS. Therefore, we recommend future studies in cycads further evaluate these loci in details.

We expect that by providing and highlighting these new resources to the plant research community, it will allow for development of new diagnostic markers and innovative conservation strategies in this iconic, but highly threatened taxonomic group, especially in the case of C. debaoensis.

Materials and Methods

DNA sequencing and genome assembly

Total genomic DNA was extracted from 0.1 g of frozen fresh leaves, from an individual collected in Guangxi (23°69′40″N, 106°15′83″E) in 2015 (voucher deposited at our research group herbarium, Jiang_C2) according to the manufacturer instructions with the Plant Genomic DNA Kit (Tiangen Biotech Co., Ltd). A 350-bp paired-end library was then constructed using NEBNext Ultra II DNA Library Prep Kit (Ipswich, Massachusetts, USA) and sequenced by Novogene (Beijing, China). About 1 Gb of raw data were obtained on an Illumina HiSeq2500 platform (San Diego, California, USA), with a paired-end read length of 2 × 150 bp. The raw reads were submitted to the SRA under the accession number SRR3407155.

The raw data were imported in Geneious R9 (Biomatters Ltd, Auckland, New Zealand), and a cp genome was assembled according to Hinsinger and Strijk51. Raw reads were trimmed according to their 5′ and 3′ -end quality, then a reference-guided assembly was performed, using the available cp of Cycas debaoensis (KM459003) as a reference for the mapping step. The cp genome annotation was transferred from C. debaoensis using the implemented function in Geneious R9 and their boundaries were manually checked. The circular cp genome map was generated using the Organellar Genome Draw program (OGDRAW).

Intergeneric comparisons and phylogenetic reconstruction

Following the recommended best practices for complete organellar sequencing52, we performed a phylogenetic analysis to confirm the accuracy of our reconstructed plastid and sample identification. We retrieved all the Cycadales available in GenBank (accessed 2016/02/15), representing all the ten genera in the order except Microcycas6. To this Cycadales dataset, we added the sequences of Gingko biloba (NC016986), Pinus strobus (NC026302), Araucaria heterophylla (NC026450), Taxus mairei (NC020321) and Cupressus sempervirens (NC026296) as outgroups.

Sequences were aligned using MAFFT53 with default options. We used the jModelTest54 implementation in CIPRES55, and set the substitution model accordingly. We built a maximum likelihood tree using PHYML56, with a 012310 + G + F model using four gamma categories and 1000 bootstrap replicates. The ML tree was built using all positions. In addition, to identify regions with substantial variability, the complete cp genomes of nine cycad genera were compared using mVISTA57, with C. debaoensis (KM459003) as a reference for the annotations.

Sequence divergences among cycads were estimated using the Kimura 2-parameter model58, implemented in MEGA659. Codon positions included were 1st + 2nd + 3rd + Noncoding. All positions containing gaps and missing data were excluded prior to analyses.

Interspecific comparisons

The complete cp genomes of C. debaoensis and two other species in cycas (C. revoluta and C. taitungensis) were compared using mVISTA57, as described above.

For each species, repeats (forward, palindrome, reverse and complement sequences) were identified using REPuter60 with 30 bp and sequence identity greater than 90%. Simple sequence repeats (SSRs) of C. debaoensis and the two other species were detected using MISA61 by setting the minimum number of repeats to 10, 5, 4, 3, 3 and 3 for mono-, di-, tri-, tetra-, penta- and hexa nucleotides, respectively. Sequence divergences among the Cycas species were estimated as described above.

Intraspecific genome comparison

The two complete chloroplasts of C. debaoensis (KU743927, KM459003) were aligned in Geneious R9 (Biomatters Ltd, Auckland, New Zealand) using the MAFFT algorithm53, and differences were identified using the “Find Variations/SNPs” function and checked individually. We recorded substitutions and indels separately, as well as their location in the chloroplast genome (e.g. SSRs/repeats, coding region/rRNA/tRNA/IGS). Sequence divergence extent between the two individuals was estimated as described above.

Additional Information

Accession codes: Raw reads and assembled chloroplast of C. debaoensis are available under the accession numbers SRR3407155 and KU743927, respectively.

How to cite this article: Jiang, G.-F. et al. Comparison of intraspecific, interspecific and intergeneric chloroplast diversity in Cycads. Sci. Rep. 6, 31473; doi: 10.1038/srep31473 (2016).