Complete plastid genome sequence of Primula sinensis (Primulaceae): structure comparison, sequence variation and evidence for accD transfer to nucleus

Library preparation and Illumina sequencing

Fresh leaves of P. sinensis were collected from the South China Botanical Garden, Chinese Academy of Sciences (CAS). Modified CTAB method was used to isolate whole-genome DNA (Porebski, Bailey & Baum, 1997). RNAs in the initial extracts was digested by RNase A to acquire pure genomic DNA. Eight primer pairs involved in this study were designed according to Yang, Li & Li (2014), to amplify the whole chloroplast genome sequence. Primers were designed to cover inverted repeat region only once. A region of approximately 16 kb was amplified by each primer pair. Long-range PCR was performed with 25 μl reaction system by using Primestar GXL DNA polymerase (TaKaRa Bio., Dalian, China). Reactions were initially denatured for 1 min at 95 °C, followed by 35 cycles of 10 sec at 94 °C, 1 min at 60 °C, and 15 min at 68 °C. Eventually the additional extension step performed on 5 min at 68 °C.

PCR products were mixed to build pair-end library using Nextera XT DNA Library Prep Kit (Illumina Inc., San Diego, CA, USA). PCR products mixture was fragmented into ∼300 bp size by the Nextera XT transposome. Library Sequencing acquired 2× 250 bp paired reads using Illumina Miseq Desktop Sequencer at Kunming Institute of Botany, CAS. All reads data were deposited to NCBI SRA database with an accession number SRP068226.

Plastome assembly and annotation

Reads were assembled using CLC Genomics Workbench v7.5.1 (CLC Bio., Aarhus, Denmark) after removing adaptors and trimmed low quality reads. Assembly was conducted twice separately with two different k-mer value 60 and 64. Contigs generated by assembling were subjected to BLAST searches against the complete chloroplast sequence of P. poissonii (NC_024543). Then, relative position and direction of each hitting contig were determined. Subsequently, hitting contigs were assembled manually to acquire complete chloroplast sequence in Geneious R7 (Biomatters, Auckland, New Zealand). The resulting plastome sequence was used as the reference, which was subsequently verified by remapping initial reads. Regions of four SC-IR junctions were identified by Sanger sequencing using four pair primers. Primer sequences used in this study can be found (Table S1).

Plastome annotation was performed using DOGMA (http://dogma.ccbb.utexas.edu/) (Wyman, Jansen & Boore, 2004) and compared with other Primulaceae species in alignment matrix. Annotations of each gene was adjusted to appropriate start and stop codons in accordance with the genetic codon for plant plastid. Incomplete genes identified from P. sinensis were verified by Sanger sequencing. The names of a few genes were updated according to the latest study, including ycf3 to pafI and ycf4 to pafII (Wicke et al., 2011). The annotated chloroplast genome sequence was submitted to GenBank (accession number: KU321892). Finally, a circular genome map was illustrated using OGDRAW (http://ogdraw.mpimp-golm.mpg.de/) (Lohse et al., 2013).

Sequence variance and SSR analysis

To further identify highly variable regions in chloroplast genome of Primula, sequences of P. sinensis and P. poissonii were compared according to variability (Pi). Alignment was performed with MAFFT version 7 (Katoh & Standley, 2013). The genetic diversity (Pi) were calculated in each split regions (400 bp) of alignment using DnaSP version 5 (Librado & Rozas, 2009). The adjacent split regions were overlapped each other with 300 bp. Gaps in the alignment were excluded from analysis. With the variance of recommended barcode rbcL + matK + ITS as reference (Yan et al., 2015), the top 40 regions with most variable sites, whose aligned lengths are longer than 200 bp were extracted for the next informative character (PICs) analysis. We followed the method of Shaw et al. (2005) to count manually the numbers of nucleotide substitutions and indels for each regions and plot them in Fig. 3.

To detect and locate simple sequence repeats (SSRs), GMATo v1.2 (Wang, Lu & Luo, 2013) was used to screen the chloroplast sequence of two Primula species. The parameter settings of mononucleotide and dinucleotide to hexanucleotide were at least eight repeat units and four repeat units, respectively.

RNAseq and accD gene characterizing

Total RNA was extracted using a modified CTAB method. Quantified total RNA (concentration ≥ 100 ng/μL; rRNA ratio ≥ 1.5) was delivered to Majorbio (Shanghai, China), where cDNA sequencing was performed with Illumina Hiseq4000 platform. Raw data were filtered and deposited in the Sequence Reads Archive (SRA) database under accession number SRX1665905. The cleaned reads were assembled de novo using Trinity with the default parameters to obtain 48,887 unisequences. Eight different software, namely TargetP (Emanuelsson et al., 2000), Protein Prowler (Hawkins & Bodén, 2006), BacelLo (Pierleoni et al., 2006), CELLO2GO (Yu et al., 2014), Euk-mPLoc2 (Chou & Shen, 2010), EuLoc (Chang et al., 2013), HybridGO-Loc (Wan, Mak & Kung, 2014), and Predotar (Small et al., 2004), were used to identify the subcellular location signals from N-terminal sequence of putative nuclear-encoded protein.

Phylogenetic analysis

Plastome of P. sinensis together with eight other plastomes published previously from different genera in Ericales was involved in this phylogenetic analysis, and Agrostemma githago from Caryophyllales was used as outgroup. All plastomes used in this study are available in GenBank (Table S2). A total of 57 plastid protein-coding genes were concatenated to generate three data sets according to the different strategies, such as all CDS, codon1+2 and condon3, respectively. On the other hand, the fourth datasets were generated including 30 pt-accD genes from six families in Ericales and putative n-accD gene in P. sinensis. All sequences from four datasets were aligned using the default option implemented in MAFFT version 7 (Katoh & Standley, 2013). Maximum likelihood (ML) trees was constructed with RAxML (RAxML-VI-HPC, http://www.trex.uqam.ca/) using GTR + Γ nucleotide substitution model (Stamatakis, 2006). Branch supports were assessed with 500 bp replicates.

Chloroplast genome assembly

The P. sinensis chloroplast genome was sequenced using Illumina Miseq sequencer, producing total number of 2× 250 bp pair-end reads 2,640,572. The mean coverage depth is about 3838.8×, ranging from 40× minimum to 25814× maximum. Two assemblies with different k-mer values successfully generated identically complete sequence with no gaps. Six contigs from assembly with k-mer value 60 were matched to reference plastome sequence, which was used to determine relative position and direction respectively. A new draft chloroplast genome was generated by identifying overlap regions manually. The draft genome was then checked and corrected according to quality and coverage of each base position by reads remapping (Fig. S1). The annotated genome was deposited into GenBank under the accession number KU321892.

Features of the P. sinensis chloroplast genome

The chloroplast genome of P. sinensis has a total length of 150,859 bp. It is divided into three parts: 82,064 bp of large single-copy (LSC) region, 17,725 bp of small single-copy (SSC) region, and two inverted repeat (IR) regions with 25,535 bp of one copy in length. The nucleotide composition of this genome has a GC content of 37.2%. Comparative analysis revealed that the genome structure of P. sinensis shared a high similarity structure to other Primulaceae species (P. poissonii NC_024543, Lysimachia coreana NC_026197, Ardisia polysticta KC465962) (Ku, Hu & Kuo, 2013; Son & Park, 2014) (Table 1).

A total of 112 genes were detected in this chloroplast genome, which could be classified into 78 protein-coding genes, 30 tRNA genes and four rRNA genes. Among them, seven coding genes, seven tRNA genes and four rRNA genes have two copies due to their locations in the IR regions. As we expected, split genes also exist in this plastome. Among them, 15 genes contain an intron and two genes have two introns (Table 2; Fig. 1).

In particular, the rps12 gene is interrupted into three pieces which resulted in trans-splicing because of its first exon located in LSC, while its second and third exons located in IRs (Hildebrand et al., 1988). Notably, the phenomena that two coding regions sharing overlapped sequence with different reading frames was found in psbD and psbC, atpE and atpB, and rps3 and rpl22. In addition, three genes, rps19, ndhF and ycf1, cross the LSC-IRb, IRb-SSC, SSC-IRa boundary, respectively. Furthermore, the ndhF 3′-terminal sequence shares the region in the IRb with the rest of ycf1 5′-terminal sequence, while the IRb-SSC boundary of P. poissonii was separated from the start codon of ndhF with 10 bp length (Fig. 2). Significantly, the accD in P. sinensis were identified as pseudogene with an extremely reduced ORF. Meanwhile, the infA lacked intact ORF in P. sinensis, P. poissonii and L. coreana, which strongly indicated that pseudogenization remains occurred in certain angiosperm groups. Further studies are necessary to focus on the mechanism of occurrence of these pseudogenes and applications of orthologous genes for phylogenetic analysis.

SSR analysis of Primula sinensis

Perfect SSRs were screened in P. sinensis and P. poissonii conducted by GMATo v1.2. Mono-, di- and tri-nucleotide repeats were found in both species (Table 3). The total number of SSRs in P. sinensis chloroplast genome is 193, of which 148 homopolymers, 44 dipolymers and 1 tripolymers are respectively found in this genome. Among them, 85.49% SSRs are only composed of A/T bases. Similar quantity level and base proportion of SSRs were also found in P. poissonii (Tables S3–S5).

Sequence variation in two Primula plastomes

Although the chloroplast genome of P. sinensis has a similar genome structure in gene contents and orders in comparison with P. poissonii, there are considerable differences in noncoding regions, especially in intergenic sequence (IGS) regions. Highly divergence regions are potential molecular genetic markers for population genetics studies. We therefore compared the regional divergence of chloroplast genome sequences of these two Primula species. The Pi value generated by DnaSP version 5 was used to indicate the level of divergence between P. sinensis and P. poissonii (Fig. S2). Then, the top 40 regions with most variable sites, with aligned length longer than 200 bp, were extracted from the alignment for further analysis (PIC calculation).

The result shows the genetic diversity (Pi) varied from 0–0.47 (Table S6). Most of the variation occurs in the non-coding regions of the LSC and SSC, while less variable characters were found in IRs. The four most variable loci, namely rpl36–rps8, rps16–trnQ, trnH–psbA and ndhC–trnV, were identified. All of these loci have been reported previously (Shaw et al., 2007). Remarkably, rpl36–rps8 exhibits the improved degree of variation. A total of 41 regions had been extracted to calculate their PICs (Fig. 3). Comparing the result with sequence variance, PICs are considered to be affected by sequence length apparently. For example, the trnH-psbA spacer has the lower number of PICs as its aligned length is only 213 bp. In contrast, ycf1 have the high PICs with the sequence length coming up to 2,500 bp in total. However, difficulties in primer designing and PCR amplification limit them (such as ycf1) for further phylogenetic use. We therefore recommend rpl32–trnL, trnS–trnG, pafII–cemA, trnC–petN, trnT–trnL, trnK–rps16 as efficient chloroplast markers, considering their balanced PICs and length size, although additional researches of their utilities are necessary in future.

Pseudogenization of pt-accD and evidence for transfer to nucleus

AccD gene encodes β subunit of the Acetyl-CoA Carboxylase, which is unique related to fatty acid synthesis within the chloroplast. Gene knockout experiments in tobacco indicated that accD gene function is indispensable, suggesting it should be an essential gene (Kode et al., 2005). So far, several researches have reported that lack of a plastid accD (pt-accD) gene or pseudogenization of pt-accD are widely present in multiple distant lineages, including Acoraceae, Campanulaceae, Ericaceae, Fabaceae, Geraniaceae and Poaceae (Goremykin et al., 2005; Guisinger et al., 2008; Haberle et al., 2008; Magee et al., 2010; Fajardo et al., 2013; Harris et al., 2013; Martínez-Alberola et al., 2013), which implied the lack of pt-accD gene or events of pseudogenization occur independently. Owing to its fundamental function in plastid development, there should be some equivalent genes replaced functionally in other subcellular structures. In recent years, a few case of nuclear encoded accD genes (n-accD) originated from plastid have been found in different taxa, such as Sciadopitys verticillata (Sciadopityaceae) (Li et al., 2016), Trifolium repens (Fabaceae) (Magee et al., 2010), Trachelium caeruleum (Campanulaceae) (Rousseau-Gueutin et al., 2013). These n-accD genes share a similar 3′-terminal region with pt-accD genes which corresponds to the carboxylase domain, while the 5′-terminal regions are completely different from pt-accD. Products of these n-accD have a putative transit peptide at the 5′-terminus. The functional prediction suggests that the transit peptide guides the protein product relocation to chloroplast. Fluorescence microscopy showed that n-ACCD-GFP fusion protein was imported in plastids contained in tobacco guard cells (Rousseau-Gueutin et al., 2013), which provided strong evidence that nuclear encoded accD gene still returns to plastid to play its roles as the same as other subunits of nuclear-encoded plastid ACCase do.

In this study, plastid accD locus in Primula sinensis shows the truncated gene with incomplete ORF. We found a portion of pt-accD gene near the start position was absent with about 400 bp in length by comparing with other plastomes from Primulaceae. This deletion has been verified by PCR amplification using a pair of primers located on both flanks. Reads mapping also showed high coverage levels on flanking regions of the deletion. Due to the presence of the deletion, pt-accD coding sequence was terminated prematurely with the introduction of the stop codon. The remaining ORF with residual sequences does not include the conserved functional region. Actually, pt-accD has become a pseudogene, with a deletion involved (Figs. S3 and S4).

Owing to its important function in plastid development, accD should probably be transferred into other subcellular structures, and still retained its catalytic activity (Rousseau-Gueutin et al., 2013). We therefore checked the transcriptome of P. sinensis using P. poissonii pt-accD as reference. The result of blast searching showed only one transcript, which was highly identical to the plastid accD functional region of P. poissonii. This transcript contained an intact ORF encoding a protein with 362 amino acids. In comparison with pt-accD sequences of other related species (such as Androsace bulleyana, Lysimachia coreana and Ardisia polystica, Fig. 4), the ORF of this transcript is different from pt-accD coding regions of other species. However, their C-terminal regions are conserved, which all contained three putative motifs, such as acetyl-CoA binding sites, carboxybiotin binding sites and carboxyltransferase catalytic sites (Lee et al., 2004). In contrast, N-terminal of the ORF shows no similarity against any plastid-encoded sequences (Fig. 4). In order to test whether the software (TargetP and Protein Prowler) affect the results of subcellular localization prediction of this putative protein-encoded sequence, we also predicted the protein sorting signals using other six software (BacelLo, CELLO2GO, Euk-mPLoc2, EuLoc, HybridGO-Loc, and Predotar). Our conclusion of the protein subcellular localization was confirmed by identical results provided by other software (data not shown). Prediction results showed a chloroplast transit peptide of 72 residues length located at the N-terminus (Fig. 4). The chloroplast transit peptides of nuclear-encoded plastid proteins (NUPTs) are necessary for targeting and import of proteins into chloroplasts. It is strongly suggested that plastid accD has been functionally transferred to the nucleus in P. sinensis. This is the fourth report in angiosperms (to our knowledge) for the transferability of accD gene with lines of solid evidences.

Phylogenetic analysis of the n-accD transcript in P. sinensis with pt-accD from other 30 species, which belongs to six families in Ericales, was performed by RAxML. This dataset contained multiple alignment of 31 sequences of the C-terminal functional regions, since the regions are relatively conserved as discussed above. The maximum likelihood tree showed n-accD from P. sinensis was located within the clade of the Primulaceae with P. poissonii as the sister group (Fig. 5). Considering the close relationship between P. sinensis and P. poissonii, this result probably indicated that accD of P. sinensis might have transferred to nucleus recently. Available data showed that complete pt-accD genes have been lost functionally from three species, Arbutus unedo (JQ067650), Vaccinium macrocarpon (JQ757046), and Chamaedaphne calyculata (KJ463365) in Ericaceae. It is likely that the missing or pseudogenization of pt-accD genes occurred accidently and independently in Ericales.

Phylogenetic analysis based on plastome sequence

The family interrelationships within the large order Ericales, constrained by the data available, have remained unclear (Anderberg, Rydin & Källersjö, 2002; Bremer et al., 2002; Geuten et al., 2004). Recently, advances in high-throughput sequencing have provided a large amount of data, which was an improvement of the phylogenetic resolution (Wen et al., 2013; Yang et al., 2015). With the expectation, nine plastome sequences represented different genera in Ericales involved in phylogenetic analysis, with Agrostemma githago as outgroup. Phylogenetic relationships were inferred using 57 plastid protein-coding genes. As we know, the third base position of codon evolves faster than the rest of two positons with higher substitution rate. ML trees were produced by using three different datasets, all coding sequence, codon 1 + 2 and codon3, respectively. The result showed that topologies of three datasets were highly congruent with one another and all nodes were well supported (Fig. 6). Primuloideae and Myrsinoideae within Primulaceae s. l. fall in two branches separately with strong supports. Primulaceae s. l. was placed as sister to the clade that comprising Theaceae, Actinidiaceae, and Ericaceae. The phylogenetic positions of these groups are in agreement with recent studies (Stevens, 2012; Zhang et al., 2015).