1. Introduction

Burmanniaceae are currently placed in the monocot order Dioscoreales, which molecular data indicates consists of Nartheciaceae, Burmanniaceae, and Dioscoreaceae (Merckx et al., 2006). However, the phylogenetic relationships between these families remains poorly understood. Burmannia is the largest genus within the family Burmanniacea, and consists of about 60 herbaceous species which use a range of trophic strategies, from autotrophy, hemi-mycohetertrophy to holo-mycoheterotrophy (Jonker, 1938; Merckx et al., 2006; Bolin et al., 2017). The genus is mainly distributed in pantropical regions, although some species can extend into warm temperate regions (Zhang and Saunders, 2000). Burmannia disticha L., is an annual chlorophyllous herbs that grows up to 12.5-70 cm tall. Stems of B. disticha are green and cauline leaves are linear. The plant commonly grows in wet thickets or grasslands and is distributed mainly in southern and southeastern Asia, and southern and south-western China (Zhang, 2001).

Whole plastid genome sequencing of heterotrophs has become increasingly common. Since the first plastid genome of a nonphotosynthetic plant was sequenced in 1992 (Epifagus virginiana; Wolfe et al., 1992), researchers have increasingly focused on plastid evolution. For instance, researchers have attempted to determine which plastid genes have been conserved in mycoheterotrophic land plant taxa by examining the plastomes of five families of monocot angiosperms (Petrosaviaceae, Orchidaceae, Corsiaceae, Triuridaceae, Thismiaceae) (Delannoy et al., 2011; Logacheva et al., 2011, 2014; Barrett and Davis, 2012; Barrett et al., 2014; Mennes et al., 2015; Schelkunov et al., 2015; Bodin et al., 2016; Lim et al., 2016).

Burmanniaceae is a small mycoheterotrophic monocot family with nearly 130 species described and in which loss of chlorophyll has occurred on at least eight independent occasions (Merckx et al., 2008). Thus, to fully clarify the evolutionary pathways from autotrophy to full mycoheterotrophy, data representing other mycoheterotrophic monocot families regarding the loss of photosynthetic activity are needed. Furthermore, previous studies have indicated that more Burmanniaceae DNA sequences are needed to clarify the phylogenetic placement of the family within Dioscoreales. Plastid phylogenomics has recently become an effective way to elucidate the difficult evolutionary relationships in plants (Mennes et al., 2013, 2015; Barrett et al., 2014; Logacheva et al., 2014; Schelkunov et al., 2015; Bodin et al., 2016; Lim et al., 2016) due to its moderate nucleotide substitution rate and suitable genome size.

In this study, we use Illumina technology and genome skimming sequencing data to characterize the complete plastome of the normal photosynthetically plant B. disticha. Specifically, we describe B. disticha plastid genome size and gene content. To explore areas of the plastid genome that might be informative in resolving Burmannia phylogeny, we compare the B. disticha genome to seven Dioscoreales genomes. This study reports the first plastid genome sequence from Burmanniaceae, and provides a reference plastid genome for plants in Dioscoreales outside of Dioscoreaceae and Nartheciaceae. Finally, this study is the first step in a more extensive comparison of plastid genome diversity within Burmannia.

2. Materials and methods 2.1. Plant material, DNA extraction, and library preparation

Plant material from B. disticha was collected in the field at Xiajinchang, Malipo County, Yunnan Province, China on October 28, 2013. The leaves were dried in silica gel and stored at room temperature for DNA extraction. Approximately 100 mg of leaf tissue was used for total genomic DNA extraction using a modified cetyltrimethyl ammonium bromide (CTAB) method (Doyle and Doyle, 1987; Rai et al., 2003). The quality of the obtained DNA was assessed by 1.0% agarose gel electrophoresis. Illumina shortinsert, paired-end (PE) libraries with average insert size of 500 bp were constructed for B. disticha, generating approximately 700 Mb clean reads of 90 bp on an Illumina HiSeq platform.

2.2. De novo contig assembly, plastome annotation, and plastome reconstruction

The genomic reads of B. disticha were assembled de novo using the SOAPdenovoV.2.0 (Luo et al., 2012) with a variety of kmer length (from 41 to 81). The resulting contigs/scaffolds were mapped to seven plastid genomes of Dioscoreales to identify potential plastid origin fragments. The process was accomplished by blasting plastid gene sequences of Aletris fauriei (GenBank accession NC_033412.1), Aletris spicata (NC_033411.1), Dioscorea elephantipes (NC_009601.1), Dioscorea rotundata (KJ490011.1), Dioscorea villosa (NC_034686.1), Dioscorea zingiberensis (NC_027090.1) and Metanarthecium luteoviride (NC_029214.1) against the assembled contig/ scaffolds pool using BLAST 2.2.26 (Altschul et al., 1990) with default settings. The aligned contig/scaffolds larger than 1000 bp were further connected by overlapped sequences, guided by their mapped locations in the reference plastid genomes. Finally, we mapped the raw sequencing reads to the assembled plastid genome sequences by Bowtie 2 (Langmead and Salzberg, 2012) to validate the assembly. The consensus sequences were annotated using the online tool DOGMA (Wyman et al., 2004) at the default settings. To complete annotation and check the gene predictions suggested by DOGMA, two additional methods were used: (1) aligning the Dioscoreales genes (as described above) using Geneious tool "MAFFT Alignment" [Version 10.1.2, Biomatters Limited; Kearse et al., 2012] and (2) annotating plastomes with the online tool CpGAVAS (Liu et al., 2012). The annotation results obtained were manually checked and corrected as necessary. In addition, we used tRNAscanSE (Schattner et al., 2005) to form more precise candidate transfer RNA (tRNA) genes. The circular plastome maps were generated in the online tool OGDRAW (Lohse et al., 2013). Finally, we determined the exact position of IR copies in sequences by aligning plastome sequences to themselves using BLASTN. A blast match was considered to describe IR if its length was more than 500 bp with identity more than 95%, and two sequences constituting the match were reverse complement to each other.

2.3. Phylogenetic trees reconstruction

We used 78 protein-coding genes from 8 complete plastomes of Dioscoreales to reconstruct the phylogenetic tree (Carludovica palmata as outgroup). C. palmata, belonging to the monocot order Pandanales, is a close, photosynthetic relative of Dioscoreales (Lam et al., 2015, 2016). Protein-coding genes were aligned by ClustalW in MEGA v 4.0 (Tamura et al., 2007), and the phylogenetic tree for a concatenated alignment was calculated in RAxML v 8.2 (Stamatakis, 2006) with a GTR+G model for the rapid Bootstrap (BS) algorithm, that was combined with the search for the best scoring maximumlikelihood (ML) tree (1000 bootstrap replicates).

3. Results and discussion 3.1. Size and structure of the B. disticha plastid genome

Illumina sequencing produced 7, 559, 832 paired-end reads, among which 698, 426 reads (9.24% out of total) could be mapped to the assembled plastid genome of B. disticha with an average sequencing depth of 399.1×. The plastid genome was assembled as a circular sequence of 157, 480 bp (GenBank accession number MG792012; Fig. 1). B. disticha has the typical quadripartite structure of plant plastid genomes, with an 81, 231-bp large single copy (LSC) region, a 13, 017-bp small single copy (SSC) region, two inverted repeat (IR) regions of 31, 616 bp, and exemplifies a typical angiosperm plastome arrangement. The overall GC content of B. disticha plastome is 34.9%. However, the distribution of nucleotides is uneven across the genome, and the IRs are more GC rich (39.5%) than the LSC (32.3%) and SSC (28.8%) regions. Approximately 120-130 unique genes are retained in the plastomes of photosynthetic land plants, and the plastid genome size is commonly about 120-170 kb (Ruhlman and Jansen, 2014). The plastome of B. disticha contains 111 genes including 78 protein-coding genes, 4 ribosomal RNA (rRNA) genes, and 29 tRNA genes (Fig. 1; Table 1). Of these, 12 genes are located in the small single copy region and 20 are duplicated in the IR. Coding regions comprise 61.3% of the B. disticha plastid genome. The mean length of intergenic regions and the longest intergenic region are 398 bp and 1918 bp, respectively, and the ratio of genic:intergenic DNA is 1.58:1.

3.2. Comparisons among the plastid genomes of B. disticha and other Dioscoreales species

Seven Dioscoreales species have plastid genomes smaller than that of B. disticha (Table 2). D. elephantipes has the smallest plastid genome (152, 609 bp) and A. spicata has the largest (154, 999 bp). The length of their LSC region varies from 81, 231 bp (B. disticha) to 85, 601 bp (D. rotundata), that of the SSC region from 13, 017 bp (B. disticha) to 19, 038 bp (D. zingiberensis), and their IR region from 25, 476 bp (D. rotundata) to 31, 616 bp (B. disticha) (Table 2). The overall GC content ranges from 37.2% to 37.5% among them, however it has been reduced to 34.9% in B. disticha. The GC content in the LSC and SSC of B. disticha are 32.3% and 28.8%, which is lower than those of the seven related plastomes, including Aletris (35.5% and 31.3%, respectively), Dioscorea (35.1% and 31.1%, respectively), and M. luteoviride (35.2% and 31.2%, respectively). With respect to the GC content in the IR, Dioscorea (43.0%) is slightly higher than M. luteoviride (42.9%) and Aletris (42.85%), but lower in B. disticha (39.5%).

The length of the IR has increased approximately 18.6-24.1% in B. disticha compared with seven Dioscoreales species. The expansion of the IR at the 3'LSC boundary to include trnH_GUG is typical for monocots and also occurs in some early diverging angiosperms (Wang et al., 2008). However, the borders in these species are not the same. Both A. fauriei and D. zingiberensis have canonical IRs that contain trnH_GUG and parts of the 5'-end of ycf1 (D. zingiberensis; 286 bp) (Fig. 2). The LSC/IRb junction (J_LB) of the remaining six Dioscoreales species (IR-expanded) is located in rps19, making 2-279 bp of the 5'-end of this gene extend into IR. By contrast, a 540-bp IR expansion into LSC is found in B. disticha, making its J_LB occur within rps19, the rps19-rpl22 intergenic region and 93 bp of the 5'-end of rpl22 extend into the IR. Lam et al. (2015) also demonstrated a further expansion of the IR at the 3'LSC boundary to include the rps19 gene in the plastid genome of C. palmata. The IRa/LSC border in all plastomes except that of B. disticha is located downstream from the trnH_GUG gene, by convention located where plastome numbering begins. However, in B. disticha the border is located downstream from the rps19 gene, placing the rps19 gene in the IRa region. Additionally, in seven related Dioscoreales species the SSC/ IR border is located in the ycf1 gene with 286-1675 bp duplicated in the IR regions. In B. disticha, the IR/SSC junction occurs within the rps15-ycf1 intergenic region, resulting in the duplication of the whole ycf1 gene sequence in IRs.

The plastome of Aletris has the most number of potentially functional genes with a total of 112, including 79 protein-coding genes, 4 ribosomal RNA (rRNA) genes, and 29 tRNA genes, whereas other Dioscoreales species, including B. disticha, contain 111 genes (Table 2). The majority of the variation in plastome genes involves the ribosome protein gene rps16. All species of Dioscoreales have lost the rps16 gene that is present in Aletris. In Aletris, a functional rps16 gene (246 bp) is located in the 2281-bp-long intergenic region between trnK_UUU and trnQ_UUG. Furthermore, the plastomes of seven related species lack a functional ycf1 gene, which has instead become a pseudogenewith onlyabout 1000 bp. In B. disticha, the two ycf1 genes in the duplicated IR region have intact reading frames and appear to be functional. Our overall results indicate that the plastid gene content of B. disticha is very similar to that of those closely related species in the order Dioscoreales.

While most start codons conformed to the commonly used canonical ATG, several atypical start codons were detected in the plastome of B. disticha, such as ACG for rpl2 and GTG for rps19. The ACG for rpl2 is presumably converted to the canonical ATG codon by RNA editing, a feature commonly found in all monocots, whether photosynthetic or not (Schelkunov et al., 2015).

3.3. Phylogenetic analyses

The final concatenated data matrix of 78 protein-coding genes was 68, 911 bp with 55, 060 (79.9%) identical sites. Phylogenetic analysis showed that B. disticha was strongly supported as the sister group to a clade that included D. elephantipes, D. rotundata, D. villosa, and D. zingiberensis (Fig. 3). Despite the potential drawbacks of phylogenies derived from plastid genomes (see Davis et al., 2014), our results are consistent with previous studies based on both mitochondrial and nuclear data which showed that Burmanniaceae is close to Dioscoreaceae, and therefore a member of Dioscoreales (Angiosperm Phylogeny Group, 1998, Angiosperm Phylogeny Group, 2003; Caddick et al., 2002a, b; Soltis et al., 2000; Davis et al., 2004; Merckx et al., 2006). Well-supported relationships in Dioscoreales include a small clade comprising M. luteoviride, A. fauriei, and A. spicata, a sister-group relationship between D. elephantipes and D. rotundata, and between D. villosa and D. zingiberensis (Fig. 3). In this study, we were unable to retrieve sufficient representative Burmanniaceae data to compare our phylogeny to those of previous studies on this family. While further taxon sampling is necessary in the future, our study provides evidence that plastome sequences may offer a new approach to resolving the phylogenetic position of Burmannia.

Acknowledgements

We thank Dr. An-Dan Zhu (Kunming Institute of Botany) for the circular drawing editing and valuable comments; This study was supported by the Large-scale Scientific Facilities of the Chinese Academy of Sciences (Grant No: 2017-LSFGBOWS-02) and the Youth Innovation Promotion Association of Chinese Academy of Sciences (2015321).