The complete chloroplast genome sequence of strawberry (Fragaria  × ananassa Duch.) and comparison with related species of Rosaceae

Plant material, DNA sequencing and genome assembly

Approximately 100 g of fresh young leaves of F. × ananassa ‘Benihoppe’ was collected from the Zhenjiang Institute of Agricultural Sciences in a Hilly Area of Jiangsu, Jurong, China. The voucher specimens were deposited in the laboratory of Fruit Tree Biotechnology of Nanjing Agricultural University. Chloroplast DNA was extracted using the high-salt saline plus Percoll gradient method of Vieira et al. (2014). A paired-end library was constructed from 50 ng of purified cpDNA according to the manufacturer’s instructions (Illumina, San Diego, CA, USA). The library, which contained an insert size of 350 bp, was sequenced using the Illumina HiSeq 2500-PE150 platform by Beijing Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). MITObim v1.8 (Hahn, Bachmann & Chevreux, 2013) was utilized for de novo genome assembly, and the chloroplast genome reads were aligned to closely related cpDNA sequences from F. vesca ssp. vesca Hawaii 4 (JF345175). Different k-mer sizes were tested, among which 31 bp produced the best results and was used to generate the final assembly in terms of the single longest scaffold length. The junctions between SC and IR regions were verified through polymerase chain reaction (PCR) amplification using sequence-specific primers (File S1). The PCR products were sequenced via Sanger sequencing.

Genome annotation and codon usage

The Dual Organellar GenoMe Annotator (DOGMA; http://dogma.ccbb.utexas.edu/, Wyman, Jansen & Boore, 2004) was employed to annotate the F. × ananassa ‘Benihoppe’ chloroplast genome. The initial annotations and putative start, stop, and intron positions were checked manually based on comparison with homologous genes in other Fragaria chloroplast genomes available in the GenBank database. Additionally, tRNA genes were identified using tRNAscan-SE 1.21 (http://lowelab.ucsc.edu/tRNAscan-SE/; Schattner, Brooks & Lowe, 2005) and ARAGORN (Laslett & Canback, 2004). A circular chloroplast genome map of F. × ananassa ‘Benihoppe’ was constructed using the online tool OGDRAW (http://ogdraw.mpimp-golm.mpg.de; Lohse, Drechsel & Bock, 2007). GC content, codon usage and relative synonymous codon usage (RSCU) were analyzed with MEGA 7 software (Kumar et al., 2008).

Repeat structure and simple sequence repeats (SSRs)

The sizes and locations of forward, reverse, palindromic and complementary repeats were determined with the REPuter program (Kurtz et al., 2001). The minimum identity and size of the repeats were limited to 90% (Hamming distance of 3) and 30 bp, respectively. SSRs in the chloroplast genome were detected using MISA (Thiel et al., 2003) with the following parameters: minimum SSR motif length of 10 bp and repeat lengths of mono-10, di-5, tri-4, tetra-3, penta-3 and hexa-3.

Comparison with other Rosaceae chloroplast genomes

One species was selected from each of the four most important fruit tree or ornamental species (Malus Mill., Pyrus L., Prunus L., and Rosa L.) of Rosaceae and from the genus Fragaria, including F. chiloensis (GP33), F. virginiana (O477), F. vesca ssp. vesca (Hawaii 4), F. vesca ssp. bracteata (MRD30) and F. pentaphylla (KY434061). The complete chloroplast genome of F. × ananassa ‘Benihoppe’ was employed as a reference and was compared with the chloroplast genomes of the nine other species using mVISTA software in the Shuffle-LAGAN mode (Frazer et al., 2004).

Nucleotide substitution in coding regions

All 78 functional protein-coding genes were extracted from the six Fragaria species (Rosoideae), Rosa roxburghii (Rosoideae) (KX768420) and Prunus persica ‘Nemared’ (Amygdaloideae) (HQ336405), and 77 protein-coding genes from Malus prunifolia (MPRUN20160302) (Amygdaloideae) (KU851961) and Pryus pyrifolia ‘Hosui’ (Amygdaloideae) (AP012207) were employed because the psbL gene was not annotated. Each exon was aligned with those of F. × ananassa ‘Benihoppe’ using ClustalX v2.1 (Thompson et al., 1997). The alignment file was then analyzed with DnaSP v5 (Librado & Rozas, 2009) to calculate synonymous (Ks) and nonsynonymous (Ka) substitution rates.

cpDNA marker identification in Fragaria

The seven complete chloroplast genomes of Fragaria species that have been released by NCBI (as above) as well as two nearly complete chloroplast genomes which were F. mandshurica (fc199s6) (KC507760) and F. iinumae (fc199s5) (KC507759) and the chloroplast genome of cultivar ‘Benihoppe’ were used to identify rapidly evolving molecular markers that may be employed for phylogenetic analysis of Fragaria. As the coding regions are highly conserved, only fragments from non-coding regions were considered. Homologous regions were aligned using MEGA 7 and adjusted manually where necessary. Then, the percentage of variable sites for each region was calculated. The proportion of mutation events = (NS/L) × 100, where NS = number of nucleotide substitutions, and L = aligned sequence length (Li et al., 2013). Because parsimony-informative sites (PIS) are commonly used in phylogenetic analyses, the number of PIS was calculated as well.

To examine the phylogenetic applications of rapidly evolving molecular markers, the maximum parsimony (MP) method was employed to construct phylogenetic trees using MEGA 7 with the following parameters: gaps in the alignment treated as missing, 1,000 replicates for bootstrap support, and tree bisection-reconnection (TBR) branch swapping.

Chloroplast genome assembly, organization, and gene content

In total, 276 Mb of 150-bp raw paired-end reads was retrieved and trimmed, and 241 Mb of high-quality short reads was finally employed to assemble the chloroplast genome, using a combination of the MITObim v1.8 de novo assembly and reference-guided (GenBank: JF345175) mapping of contigs to generate complete chloroplast genome sequences. Finally, the generated data were assembled into the single longest scaffold spanning the F. × ananassa ‘Benihoppe’ chloroplast genome. To validate the assembly, four junctions between SC and IR regions were confirmed through PCR amplification and Sanger sequencing. No mismatches or InDels were observed between the Sanger sequencing and the assembled genome, which verified the correctness of our genome sequencing and assembly results.

The F. × ananassa ‘Benihoppe’ chloroplast genome is a typical circular double-stranded DNA molecule with a quadripartite structure; it is 155,549 bp in size and consists of IR (25,936 bp) regions separated by LSC (85,531 bp) and SSC (18,146 bp) regions (Fig. 1, Table 1). The GC content of the chloroplast genome is 37.23% (Table 1), and the GC contents of the LSC and SSC regions (35.12% and 31.14%) are lower than those of the IR regions (42.85%). The high GC contents in the IR regions are mainly due to the high GC contents of the four ribosomal RNA (rRNA) genes (55.43%) which is similar to most of other plants cp genomes (Wang, Shi & Gao, 2013; Shen et al., 2016; Kong & Yang, 2017).

There are 130 genes in the chloroplast genome of F. × ananassa ‘Benihoppe’, 112 of which are unique, including 78 protein-coding genes, 30 tRNA genes and 4 rRNA genes. The genes that are repeated in IRs comprise seven protein-coding genes, seven tRNA genes, and four rRNA genes (Fig. 1, Table 2). Among these genes, a single intron was detected in 15 genes (nine protein-coding genes and 6 tRNA genes), while two genes (ycf3 and clpP) were found exhibit two introns each (File S2). The trnK-UUU gene harbors the largest intron (2,497 bp), which contains the matK gene, whereas the intron of trnL-UAA is smallest (422 bp). Two genes with internal stop codons (ycf15 and ycf68) and one without a stop codon (infA) were annotated as pseudogenes. Absence or pseudogenization of these three genes has also been reported in other Rosaceae species, such as Eriobotrya japonica (Shen et al., 2016) and Prinsepia utilis (Wang, Shi & Gao, 2013). The rps12 gene is a trans-spliced gene with a 5′ exon located in an LSC region and two 3′ exons located in IR regions. The complete chloroplast genome with gene annotations has been deposited in the NCBI GenBank database (accession number: KY358226).

Overall, 22,709 codons encoding 78 protein-coding genes were identified in the complete chloroplast genome and classified according to codon usage (File S3), among which 2,405 (10.59%) encode leucine (the most abundant amino acid), and 252 (1.11%) encode cysteine (the least abundant amino acid). RSCU analysis showed A/T contents of 53.89%, 61.84% and 70.17% at the first, second and third codon positions, respectively. This pattern of higher A/T bias at the third codon position is common in the chloroplast genomes of land plants (Morton, 1998; Gurusamy & Seonjoo, 2016).

Repeat structure and SSR loci

A total of 39 repeat structures with a minimal length of 30 bp and minimal identity of 90% were found (File S4), including 14, 6, 2 and 17 forward, reverse, complementary, and palindromic structures, respectively. Among these structures, the longest is 67 bp and is located between trnM-CAU and atpE. Most of the repeat structures are located in intergenic regions (65.4%), while fewer than half are located in coding genes (21.8%; ndhA, ycf2, psaA, psaB, trnS-GGA, psbJ, trnG-UCC and trnG-GCC) or introns (12.8%; ycf3, ndhB and clpP).

Simple sequence repeats (SSRs) in chloroplast genomes have become valuable molecular markers because of their high degree of variation within an individual species, which is useful for linkage map construction and plant breeding (Powell et al., 1995; Xue, Wang & Zhou, 2012). In the F. × ananassa ‘Benihoppe’ chloroplast genome, 61 SSR loci with a length of at least 10 bp were detected, among which 38 (62.3%) are mononucleotide repeats; 16 (26.2%) are di-repeats; three (4.9%) are tri-repeats; and four are (6.6%) tetra-repeats. No pentanucleotides or hexanucleotides were found. Most of the observed mononucleotide repeat sequences consist of A/T motifs, whereas only one is composed of a G/C motif. Similarly, 93.75% of the dinucleotide repeat sequences consist of AT/TA motifs. The results showed that the SSRs exhibit a strong AT bias, which is consistent with other studies (Lei et al., 2016; Kuang et al., 2011). Among the 61 SSR loci, 44 are located in intergenic regions, eight in introns, and nine in coding regions of genes (Table 3).

Comparison with other chloroplast genomes from Rosaceae

Nine chloroplast genomes representing five genera in Rosaceae were compared with that of F. × ananassa ‘Benihoppe’ (Table 1). The length of the Fragaria chloroplast genomes ranges from 155,549 to 155,691 bp, with F. vesca ssp. vesca (Hawaii 4) exhibiting the largest chloroplast genome and F. × ananassa ‘Benihoppe’ the smallest. The length of the LSC regions shows greater variation, ranging from 85,531 to 85,605 bp, with F. vesca ssp. vesca (Hawaii 4) exhibiting the longest, followed by F. virginiana (O477), while F. × ananassa ‘Benihoppe’ harbors the shortest (Table 1). However, the IR regions of diploid strawberries are longer than those of three octoploid strawberries. F. pentaphylla (KY434061) exhibits the shortest SSC regions. Furthermore, the size of the F. × ananassa ‘Benihoppe’ chloroplast genome is smaller than those of the other four species in Rosaceae, being approximately 4.5 kb, 4.4 kb, 2.2 kb and 1.2 kb smaller than those of M. prunifolia (MPRUN20160302), P. pyrifolia ‘Hosui’, P. persica ‘Nemared’, and R. roxburghii (KX768420), respectively. The differences in genome size can largely be attributed to variation in the length of SSC and IR regions (Table 1).

The results also revealed that the gene content and gene order of F. × ananassa ‘Benihoppe’ are identical to those of the five previously reported the genus Fragaria chloroplast genomes. Interestingly, the loss of a group II intron of the atpF gene, as observed in Fragaria (Table 2), has previously been reported for Malpighiales (Daniell et al., 2008) and R. roxburghii (KX768420). However, the numbers of unique genes found in the F. × ananassa ‘Benihoppe’, R. roxburghii (KX768420), M. prunifolia (MPRUN20160302), P. pyrifolia ‘Hosui’, and P. persica ‘Nemared’ chloroplast genomes were 112, 114, 111, 111 and 112, respectively, due to the absence of the psbL gene in M. prunifolia (MPRUN20160302) and P. pyrifolia ‘Hosui’, the absence of the trnG-GCC gene in R. roxburghii (KX768420), and the presence of three genes, infA, trnP-GGG and trnM-CAU, only in R. roxburghii (KX768420). The GC content among the ten species was similar, ranging from 36.56 to 37.25%, with the seven Rosoideae species all exhibiting a high GC content, of approximately 37.2% (Table 1).

The mVISTA program was employed to analyze the overall sequence identity among all ten Rosaceae members at the chloroplast genome level, using the annotation for F. × ananassa ‘Benihoppe’ as a reference (Fig. 2). The results showed high similarity among the Fragaria chloroplast genome sequences, particularly for F. × ananassa ‘Benihoppe’, F. chiloensis (GP33) and F. virginiana (O477). Among the other Rosaceae species, the F. × ananassa ‘Benihoppe’ chloroplast genome was most similar to that of R. roxburghii (KX768420) and most divergent from that of P. persica ‘Nemared’. Overall, the results revealed SC regions to be more divergent than IR regions, with higher divergence being observed in non-coding regions than in coding regions, which is a common phenomenon in the chloroplast genomes of angiosperms (Yao et al., 2015; Ni et al., 2016; Asaf et al., 2016). The coding regions with marked differences include the ycf1, matK and psaI genes. The highest divergence in non-coding regions was found for rps16-trnQ, petN-psbM, ndhC-trnV, petA-psbL and rpl32-ccsA. These results are similar to those of other analyses performed in Rosaceae (Wang, Shi & Gao, 2013; Shen et al., 2016), suggesting that these regions evolve rapidly in Rosaceae.

IR contraction and expansion

In general, IR regions are considered to be the most conserved regions in the chloroplast genome. Nevertheless, expansion and contraction of the border region between SC and IR regions are common during evolution and contribute to variation in chloroplast genome length (Wang et al., 2008; Li et al., 2013). Thus, the positions of LSC/IRA/SSC/IRB borders and the adjacent genes in the ten Rosaceae chloroplast genomes were aligned (Fig. 3). The SSC/IRB boundary of F. × ananassa ‘Benihoppe’ is consistent with those of the other Fragaria species. All of the genomes except for those of F. vesca ssp. vesca (Hawaii 4) and F. pentaphylla (KY434061) exhibit IRA/SSC boundaries of the same length, and due to contraction of the IR region at the IRB/LSC boundary, F. ×  ananassa ‘Benihoppe’ exhibits the shortest IR region among the six Fragaria species.

Compared with those of other Rosaceae species, the rps19 genes of Fragaria species and R. roxburghii (KX768420) are shifted to an LSC region with a 12–21 bp gap. However, the rps19 genes of M. prunifolia (MPRUN20160302), P. pyrifolia ‘Hosui’ and P. persica ‘Nemared’ extend from the LSC to the IRA region, showing variability of 120–182 bp, resulting in the presence of an rps19 pseudogene of the same length in IRB. The SSC/IRB boundary extends to the ycf1 coding region, ranging from 1,051 bp (P. persica ‘Nemared’) to 1,106 bp (R. roxburghii, KX768420), leading to a nonfunctional ycf1 gene in IRA. The ndhF gene is located entirely in the SSC region in Rosoideae species but varies in distance from the IRA/SSC border. However, in M. prunifolia (MPRUN20160302), P. pyrifolia ‘Hosui’, and P. persica ‘Nemared’, part of the ndhF gene is located in IRA. In general, the position of the trnH gene in the chloroplast genome is quite conserved between monocot and dicot species. In monocots, the trnH gene is located in the IR region, whereas it is located in the LSC region in dicots (Asano et al., 2004). In all of the analyzed genomes, the trnH gene is located in the LSC region, although its distance from the IRB/LSC junction ranges from 0 to 101 bp. Overall, a similar pattern of expansion and contraction of IR/SC regions was observed among the Fragaria species and R. roxburghii (KX768420), differing from M. prunifolia (MPRUN20160302), P. pyrifolia ‘Hosui’ and P. persica ‘Nemared’ (Fig. 3).

Selection pressure on the F.× ananassa ‘Benihoppe’ chloroplast genome

The Ka/Ks ratio was calculated for 78 protein-coding genes in all nine chloroplast genomes, with a value of 0 indicating neutral selection. The Ka/Ks ratio of the Fragaria chloroplast genomes was typically calculated to be 0, except for six genes in F. vesca ssp. vesca (Hawaii 4) (rpoC2, ndhD, ndhF, psbB, ycf1 and ycf4), three genes in F. vesca ssp. bracteata (MRD30) (ndhF, ycf1, and ycf4) and thirteen genes in F. pentaphylla (KY434061) (rpoC1, rpoC2, atpB, atpH, ndhA, ndhD, ndhF, ndhH, petA, psbB, rbcL, ycf1 and ycf4) (File S5).

Among the protein-coding genes in the chloroplast genomes of the other Rosaceae species, the Ka/Ks ratio was observed to be highest in genes within the SSC regions (Table 4). In the comparison of Fragaria with Rosa and Malus, the lowest Ka/Ks ratio was found in the IR region. However, in the comparison of Pyrus and Prunus, the LSC region showed the lowest Ka/Ks ratio (Fig. 4, Table 4). The lowest Ka/Ks ratio was observed for genes encoding subunits of ATP synthase, subunits of the cytochrome b/f complex, subunits of photosystem II and the large subunit of RuBisCO (File S5). With the exception of the rpl16 gene of Rosa, the Ka/Ks ratio of all genes was found to be less than 1, suggesting purifying selection on these genes (Fig. 4).

Variation in chloroplast DNA in three octoploid strawberries

The complete chloroplast genomes were found to be most similar among the three octoploid strawberries. However, when the sequences of the coding and non-coding regions of F. × ananassa ‘Benihoppe’ were compared in detail with those of F. chiloensis (GP33) and F. virginiana (O477), a number of SNPs and InDels were revealed (Table 5).

In total, 35 SNPs (26 transversions and 9 transitions) were identified between the complete F. × ananassa ‘Benihoppe’ and F. chiloensis (GP33) chloroplast genomes, which were found in all types of regions (23 in LSC, 9 in SSC and 3 in IR regions). Two SNPs in the rpoB and ndhF genes represent synonymous changes, whereas the other six SNPs in four other genes (accD, ndhH, ycf1 and ycf4) are nonsynonymous and may alter the encoded protein’s primary structure. Overall, 18 InDels between 1 and 19 bp in length were found, including 15 within LSC regions. In contrast, 23 SNPs (17 transversions and 6 transitions) were identified between the F. × ananassa ‘Benihoppe’ and F. virginiana (O477) chloroplast genomes, among which 21 are located in LSC and 2 in SSC regions, while none were found in IR regions. Three nonsynonymous SNPs were found in the petB, ycf1 and ycf4 genes, while three synonymous SNPs were found in rps8, rpoB, and psbA. Eight InDels were observed (5 insertions and 3 deletions), all but one of which is located within the LSC region.

Regardless of location or subspecies, all F. virginiana individuals share the same SNPs in two genes, petD (G) and ndhF (A), differing from F. chiloensis at these positions (petD: A; ndhF: (G) (Salamone et al., 2013). Furthermore, Honjo et al. (2009) examined two non-coding regions and concluded that cultivar ‘Benihoppe’ exhibits haplotype V, which is consistent with its female parent Akihime. Based on our SNP analysis, two genes (petD and ndhF) and two non-coding regions (trnL-trnF and trnR-rrn5) are consistent with F. virginiana and different from F. chiloensis. Our results are in accordance with previous studies (Salamone et al., 2013; Honjo et al., 2009) and indicate that the sequence identity between F. × ananassa ‘Benihoppe’ and F. virginiana (O477) at the chloroplast level is higher than between F. × ananassa ‘Benihoppe’ and F. chiloensis (GP33).

cpDNA markers and sequence polymorphisms in Fragaria

Non-coding regions (introns and intergenic spacers), harboring more sequence divergence, are not subject to the functional constraints that could extend the utility of a molecule at lower taxonomic levels (Small et al., 1998; Shaw et al., 2007). At least six non-coding regions of cpDNA were previously examined in phylogenetic and ancestry studies of Fragaria. Universal primers targeting the trnL-trnF region (Taberlet et al., 1991) were used in previous studies to investigate 14 species of Fragaria (Potter, Luby & Harrison, 2000) and 8 diploid species from this genus (Sargent, 2005). The trnT-trnL, atpB-rbcL, psbA-trnH, psbJ-psbF and rps18-rpl20 regions have also been employed to detect sequence polymorphisms in Fragaria (Sargent, 2005; Lin & Davis, 2000). However, due to the small group of taxa sampled or a low level of sequence variation in these regions, the phylogenetic resolution within Fragaria has been limited (Rousseau-Gueutin et al., 2009).

To examine which regions might be applied for Fragaria phylogenetic analysis, all of the non-coding regions among ten Fragaria chloroplast genomes were aligned, and sequence divergence was calculated (Fig. 5). The results showed that ndhE-ndhG exhibits the highest rate of variation (7%), while atpF-atpH and trnT-psbD exhibit the most PIS. However, only the 6 intergenic regions (trnK-matK, trnS-trnG, atpF-atpH, trnC-petN, trnT-psbD, and trnP-psaJ) display a percentage of variable sites higher than 1% and more than five PIS (Fig. 5), indicating the low variation of chloroplast genomes in Fragaria. Interestingly, these intergenic regions are all located in the LSC region, whose sequence has been noted to be less conserved in those of IR and SSC regions and has consequently been used for phylogenetic analysis at low taxonomic levels (Njuguna, 2010).

To examine the phylogenetic applications of the six fast-evolving DNA regions, an MP tree was constructed for each molecular marker from the ten Fragaria species (File S6). The results revealed that none of each region was efficient in resolving the relationships among the examined samples. However, the combined regions strongly supported F. iinumae (fc199s5) and F. pentaphylla (KY434061) had the closest phylogenetic relationship. Furthermore, our results showed F. vesca ssp. vesca (Hawaii 4) was closer to F. vesca ssp. bracteata (LNF40), which was similar to Govindarajulu et al. (2015). The STEMhy and PhyloNet results showed a greater contribution of F. iinumae than F. vesca to the ancestry of the octoploids (Kamneva et al., 2017). Vining et al. (2017) used POLIMAPS to resolve F. × ananassa chromosomal regions derived from diploid ancestor F. vesca. Our results couldn’t infer which one was the ancestor of F. × ananassa ‘Benihoppe’ (KY358226) at present. Further studies with a broad sampling scheme need to be conducted to test the efficiency of these six identified regions in phylogenetic analysis of Fragaria.