The complete chloroplast genome sequence of an endemic monotypic genus Hagenia (Rosaceae): structural comparative analysis, gene content and microsatellite detection

DNA extraction and sequencing

Young leaf samples were collected from natural populations of Hagenia abyssinica in Mt. Kenya (Kenya; 00°09′35.29″S/037°26′56.40″E). A voucher specimen (SAJIT_001956) was deposited at the Herbaria of Wuhan Botanical Garden, Chinese Academy of Sciences (HIB). Total genomic DNA was extracted from 100–150 mg of leaves using the MagicMag Genomic DNA Micro Kit (Sangon Biotech Co., Shanghai, China) following the manufacturer’s instructions. The quality of the extracted DNA was checked by gel electrophoresis and confirmed using Qubit DNA Assay kit in Qubit 2.0 Fluorometer (Life Technologies, San Diego, CA, USA). Paired-end library was constructed using an Illumina TruSeq Library preparation kit (Illumina, San Diego, CA, USA) following the manufacturer’s protocol. Genomic DNA was sequenced using the Illumina Hiseq 2500 platform (Illimina Inc.), yielding 41.2 million 150-bp paired-end reads from a library of ∼350  bp DNA  fragment.

Genome assembly and annotation

We used a reference-guided strategy to assemble the chloroplast genome. Firstly, whole clean data were identified using BLAST (http://blast.ncbi.nlm.nih.gov/) with default parameters, by searching against the plastome sequences of Fragaria chiloensis (JN884816). The generated contigs were sorted, and the chloroplast genome reads were extracted by mapping the contigs against already available chloroplast sequences of Fragaria chiloensis (JN884816; Salamone et al., 2013) using the Basic Local Alignment Search Tool (BLAST; http://blast.ncbi.nlm.nih.gov/) with default parameters. The retained high quality reads were then assembled into non- redundant contigs using Velvet 1.2.10 (Zerbino & Birney, 2008) with K-mer length of 95–107. Five contigs whose size ranged between 1,960 and 47,845 bp were then blasted against Fragaria chiloensis and Pyrus pyrifolia (AP012207; Terakami et al., 2012). Specific primers were designed using PRIMER 5.0 (PREMIER Biosoft International, CA, USA) and used in Polymerase Chain Reaction to fill gaps between the contigs and to validate the joints between the IR/LSC and IR/SSC, based on the Sanger sequencing technique. The primer sequences used in filling the gaps and validating the IR/SC junctions are listed in File S1.

The assembled chloroplast genome was annotated using an online-based program: the Dual OrganellarGenomMe Annotator (DOGMA; http://dogma.ccbb.utexas.edu/, Wyman, Jansen & Boore, 2004) followed by manual corrections of the start, stop codons and the boundaries between the introns and exons based on homologous genes from other sequenced chloroplast genomes. Protein coding, transfer RNA (tRNA) and ribosomal RNA (rRNA) genes were also predicted in DOGMA with default parameters. The tRNA genes were further verified using tRNAscan-SE 1.23 program (http://lowelab.ucsc.edu/tRNAscan-SE/; Schattner, Brooks & Lowe, 2005). Finally, a circular gene map was constructed using the OrganellarGenomeDRAW software (OGDRAW; http://ogdraw.mpimp-golm.mpg.de). The complete chloroplast genome sequence of H. abyssinica can be found in GenBank under the accession number KX008604.

Microsatellite discovery and comparative analyses

The Perl script based Microsatellite identification tool (MiSa) (Thiel et al., 2003) was used to detect microsatellites with minimal iterations of eight repeat motifs for mononucleotides, five for dinucleotides, four for trinucleotides and three for Tetra-, Penta- and hexa-nucleotides. The location and size of the repeating sequences (forward, reverse, complementary and palindromic) were visualized in REputer (Kurtz & Schleiermacher, 1999) with minimal repeat size set at ≥15 and Hamming distance at 3.

To highlight structural differences and similarities between H. abyssinica and other already sequenced chloroplast genomes in Rosaceae family, we retrieved 20 currently available complete chloroplast genomes from the NCBI (Table 1) and conducted comparative analyses. Special attention was paid to the sizes of the entire complete genomes and inverted repeats, the location of the IR/SC junctions and arrangement of genes adjacent the IR/SC boarders.

To gain insight into the relationship of members of Rosaceae, a Maximum Likelihood (ML) phylogenetic tree was reconstructed. We used 71 protein-coding genes common in all the 21 species of Rosaceae. Two taxa; Morus indica (Moraceae) and Eleagnus macrophylla (Eleagnaceae), from the clade Rosales, were used as outgroups. All the PCGs were aligned in MUSCLE (Edgar, 2004) with default settings and appropriately edited manually. The jModelTest 2.1.7 program (Darriba et al., 2012) was used to select the best fitting substitution model based on the Akaike information criterion (Posada & Buckley, 2004). The best-fitting substitution model GTR + I + G model of all genes was used. The GTR + I + G model was used for ML analyses implemented in RAxML 8.0.20 following instructions from the manual (Stamatakis, 2014). A bootstrap analysis was performed with 1,000 replications.

Genome content and organization

The complete chloroplast genome of H. abyssinica exhibited a double- stranded circular DNA molecule, with a total length of 154,961 bp (Fig. 1). It also displayed a quadripartite structure, typical to chloroplast genomes of most terrestrial plants. The chloroplast genome possesses a pair of inverted repeats (IRa and IRb) of 25,971 bp each. The IRs are separated by a large single copy (LSC) and a small single copy (SSC) with 84,320 bp and 18,696 bp respectively (Fig. 1). The total GC content for this chloroplast genome is 37.1%, which is consistent with those from other species in Rosaceae. The chloroplast genome of H. abyssinica encodes 129 genes (excluding the ORFs and the the hypothetical genes; ycf68 and ycf15), comprising 78 unique protein—coding genes (PCGs), 30 unique tRNA and 4 rRNA genes (Table 2). In total there were 17 duplicated genes, 7 of which code for protein in the IRs including rpl2, rpl23, ycf2, ndhB, rps7, rps12, and ycf1, 6 tRNA and 4 rRNA were also among the duplicates in the IRs. The gene order in the SSC region begins with ndhF, followed by rpl32, trnL, ccsA, ndhD, PsaC, ndhE, ndhG, ndhI, ndhA, ndhH and rps15 and ends with ycf1. Six protein coding genes contained either one intron (rps16, rpl2, rpl23, rpoC1, ndhA and ndhB) or two introns (clpP). The hypothetical gene ycf3, contained two introns (Table 2). The rps12 gene is trans-spliced with the 3′ exon being duplicated in the IR, while the 5′ end is located at the LSC region.

Discovery of SSRs

Microsatellite markers are considered ideal for plant molecular studies due to their high mutation rates, multi- allelism and locus- specificity (Varshney, Graner & Sorrells, 2005; Govindaraj, Vetriventhan & Srinivasan, 2015) and thus highly informative. Recently, seventeen species-specific nuclear SSR markers have been reported for this species (Gichira et al., 2016). In a previous study, three concensus chroloplast microsatellite markers had been used to study genetic diversity of H. abyssinica (Ayele et al., 2009). Chloroplast-derived microsatellite markers have generated great impact on population genetics, plant evolutionary studies and phylogenetics (Provan, Powell & Hollingsworth, 2001). In this study, a total of 172 SSR repeat motifs were discovered (Table 3).

Mononucleotides had the highest number of repeats (88%), most of which had the A/T repeat type which is in line with the findings of a previous study that polyA and polyT repeats dominate in chloroplast microsatellites (Cai et al., 2008). A total of 5.8% represented dinucleotides while the rest were tetranucleotides, nine of the dinucleotides had the AT/TA repeat motif while AAAT/TTTA motifs dominated among the tetranucleotides. There were no trinucleotide repeats detected in H. abyssinica’s chloroplast genome. Repeat motifs are potential molecular tools for studying recombination and rearrangement in genomes (Smith, 2002). In addition to SSRs, a total of 49 repeat sequences with at least 21 bp were identified by REPuter. The repeat units had a sequence identity of ≥90% and their sizes ranged from 21 to 69 bp. The 49 repeats constituted 22 palindrome (inverted) repeats, 19 direct (forward) repeats, seven reverse repeats and one complementary repeat (Table 4). The majority of the identified repeats were located in the non-coding regions of the genome which is in line with observations made in other chloroplast genomes of angiosperms (Provan, Powell & Hollingsworth, 2001; George et al., 2015). This trend of cpSSR distribution, has been observed in other chloroplast genomes species in Rosaceae suggesting that they may be suitable for conducting population genetic diversity, phylogenetic and evolutionary studies in species under this family.

Comparative analysis and phylogenetics

The number of species from the Rosaceae family with completely sequenced cholorplast genomes is rapidly increasing. Currently, the complete chloroplast genomes of 20 species from eight genera in two sub-families of Rosaceae family have been sequenced and deposited at GenBank (http://www.ncbi.nlm.nih.gov/). Out of the 20 species, 12 belong to the Spiraeaideae sub-family while the rest fall under the Rosoideae sub-family (Potter et al., 2007; Hummer & Janick, 2009). We compared the structure of Hagenia’s chloroplast genome to those available from the eight genera. The list of the species used for comparision and their accession numbers are shown in (Table 1). Characteristically, there are four junctions in the chloroplast genomes of angiosperms, due to the presence of two identical copies of the inverted repeats. However, the loss of one inverted repeat has been reported in some flowering plants e.g., in legumes (Palmer et al., 1987b). All chloroplast genomes appeared to be structurally similar with a typical quadripartite structure of two IRs separated by a LSC and a SSC. The whole genome sizes ranged from 154,959 (Potentilla micrantha) to 160,041 (Malus prunifolia) and there was a clear distinction of the sub-families based on genome sizes. Species from the Maloideae sub- family have a larger chloroplast genome compared to those from the Rosoideae. The size of H. abyssinica’s chloroplast genome (154,961 bp) is only 2 bp larger than that of the smallest chloroplast genome of P. micrantha (154,959 bp; Ferrarini et al., 2013).

Size variations of the chloroplast genome may be attributed to the expansion/contraction of the IR, with small variations (<100 bp) being common even among species under the same genus (Goulding et al., 1996). The expansion and/or contraction of the IRs is regarded as a significant evolutionary event and can be a source of polymorphic genetic markers for species identification and for analyzing phyologenetic studies in plants (Wang et al., 2008). In this study, sizes of the IRs varied from 26,435 bp in Prunus maximowiczii to 25,530 in P. micrantha. Although certain genes near the IR/SC boarders appeared to be conserved in all the species, key variations were noted in gene arrangement along the IR/SC junctions (Table 1). Two genes (rps19 and rpl2) are adjacent the IRa/LSC boarder at varying positions, while the IRb/LSC junction is flanked between genes rpl2 and trnH-GUG and in some cases a pseudogene (Ψ) of rps19 gene is included in this region. This is a common feature in angiosperms, excluding monocots whose trnH-GUG gene is located in the IR between the genes rpl2 and rps19 (Goulding et al., 1996; Wang et al., 2008).

In all species from the Spiraeoideae subfamily and in one Rosoideae species—Pentactina rupicola—the IRa/LSC junction occurs within the coding region of the rps19 gene resulting into the presence of Ψrps19 gene of various length in the IRb. This event has also been reported in the chloroplast genomes of other species e.g., Arabidopsis thaliana (Sato et al., 1999) and Coffea arabica (Samson et al., 2007). However in the other species, including H. abyssinica the entire rps19 gene is located in the LSC region, leaving a gap of varying length between the 5′ end of the gene and the IRa/LSC junction, this is similar to other dicots such as Nicotiana tabacum (Shinozaki et al., 1986). The largest gap was 1,016 bp in P. micrantha followed by 130 bp in H. abyssinica. The rpl2 gene is entirely located in both IRs region in all species, consequently leaving a gap of non-coding region between the IR/LSC junction and rpl2 gene. The IRb/LSC junction is situated in the down-stream of non-coding region of the trnH-GUG gene in all analysed species. Those species with Ψrps19, the pseudogene was located within the IR, between the rpl2 and the trnH-GUG. In some dicots e.g., Actinidia chinensis (Yao et al., 2015), trnH-GUG and a section of the psbA occur in the inverted repeat due to expansion on the IRs into the LSC region.

In all the studied species, the IRb/SSC junction is located within the coding region of the ycf1 gene. Consequently, the ycf1 gene extends into the IRb at varying lengths ranging from 946 bp in Prunus persica to 1,091 bp in all species of genus Fragaria. As a result, the IRa/SSC junction is bordered by Ψycf1 and gene ndhF, which is a general structure among the dicots e.g., tobacco and Arabidopsis. In Hagenia, the ycf1 gene has an extension of 1,040 bp into the IRb and therefore, its Ψycf1 of 1,151 bp overlaps with ndhF (2,234 bp) at 65 bp. The chloroplast genome of Annona cherimola, which is one of the largest plastid genomes with 201,723 bp, has an extremely reduced SSC (2,966 bp) due to major expansions of the IRs and most genes including the ycf genes have been incorporated in the IRs (Blazier et al., 2016).

Chloroplast DNA is reported to have evolved from free-living Cyanobacteria through endosymbiosis with a history of more than 1.2 billion years and since then a number of genes, initially found in the chloroplast genomes have relocated to the nuclear genome (Timmis et al., 2004), e.g., in Arabidopsis 18.1% of its functional nuclear genes originated from the plastid genome (Martin et al., 2002). Further studies presented more evidence on independent gene transfers from the chloroplast to the nuclear genome in rosids (Millena et al., 2001), these includes the successful transfers of rpl22 gene in Castanea, Quercus and Passiflora (Jansen et al., 2011),  infA gene in Arabidopsis (Sato et al., 1999) and in Elaeagnus (Choi, Son & Park, 2015). These transfers occurred in the initial stages of plastid evolution, though a high relocation rate of non-coding DNA happens continuously (Martin et al., 2002; Timmis et al., 2004). Generally, loss and/or transfer of genes from the chloroplast genomes to the nuclear or mitochondria genomes is as a result of evolutionary events, allowing chloroplast genomes to act as valuable molecular tools in phylogenetic and evolutionary studies. Further comparative analyses revealed that the initiation factor 1 (infA) gene which was observed in other species of Rosaceae, is conspicuously missing from the Hagenia chloroplast genome. The loss/transfer of the infA gene, which is an essential gene in Escherichia coli (Cummings & Hershey, 1994), is common among the angiosperms and it is regarded as a highly mobile gene (Millena et al., 2001; Daniell et al., 2016). Therefore, besides the expansion/contraction of the IRs, gene loss provides crucial information that is essential for evolutionary studies and resolution of phylogenetic relationships among plant species.

Complete chloroplast genome sequences provide essential genetic data for precise systematics and phylogenetic resolutions in plants. The ML phylogenetic tree that was constructed using 71 PCGs, common in all 21 taxa from Rosaceae and in two outgroups, clearly placed the Rosaceae species into two clades. The two main clades concurred with two sub-families: Spiraeoideae and Rosoideae (Fig. 2). This classification was in agreement with the phylogeny of Rosaceae (Potter et al., 2007). Previously, Hagenia had been classified in sub- family Rosoideae under Agrimoniinae, a subtribe in the tribe Sanguisorbeae, alongside the genera Aremonia, Agrimonia, Leucosidea and Spenceria  (Eriksson et al., 2003; Potter et al., 2007).