For each species, 2 g of fresh leaves was harvested from individuals of Cypripedium japonicum Thunb, Dendrobium huoshanense, D. loddigesii Rolfe, D. moniliforme (L.) Sw., Goodyera schlechtendaliana Rchb. f., Neuwiedia zollingeri var. singapureana (Wall. ex Baker), Paphiopedilum armeniaclum, and Vanilla aphylla Blume that had been cultivated in the greenhouse of Nanjing Normal University. Total DNA was extracted using a Qiagen DNeasy plant mini Kit (Qiagen, Germany). The quality of obtained DNA was measured on a NanoDrop 8000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA). DNA samples that passed the threshold (DNA concentrations greater than 300 ng/μL, A₂₆₀/A₂₈₀ = 1.8-2.0, and A₂₆₀/A₂₃₀ larger than 1.7) were collected for next-generation sequencing.

With the development of next-generation sequencing technologies, many high-throughput methods can now be used to obtain complete plastome sequences, using the Illumina platform whole-genome sequence (e.g., Kim K. et al., 2015). In this study, plastome assembly was done using the CLC Genomics Workbench 6.0.1 (CLC Bio, Aarhus, Denmark). For each species, approximately 3.7 Gb 73 bp pair-end reads were sequenced on an Illumina Genome Analyzer 2000 at Beijing Genomics Institute. After the adaptors were removed, reads were trimmed with an error probability < 0.05 and de novo assembled on the CLC Genomics Workbench 6.0.1. Contigs that had >30× sequencing depths were searched using blastn against the plastomic sequences of G. fumata (NC_026773). Matched contigs with E values < 10^-10 were designated plastomic contigs. For D. moniliforme, G. schlechtendaliana, P. armeniaclum, and V. aphylla, gaps between plastomic contigs were closed using sequences of amplicons obtained via PCR with specific primers. The boundaries of IRs were confirmed by PCR assays. Genes were predicted using DOGMA (Wyman et al., 2004) and tRNAscan-SE 1.21 (Schattner et al., 2005). The exact boundaries of predicted genes were confirmed by aligning them with their orthologs from other orchid species.

The comparative plastomic method, which uses at least two complete plastomes within the study genus, is now available for mutational hotspot selection (e.g., Ahmed et al., 2013; Shaw et al., 2014). However, only 10 orchid genera contain two or more plastomes (2 in Bletilla, 5 in Corallorhiza, 11 in Cymbidium, 2 in Dendrobium, 2 in Masdevallia, 3 in Phalaenopsis, 4 in Goodyera, 2 in Cypripedium, 2 in Paphiopedilum, and 2 in Vanilla) (Supplementary Table S1). Therefore, sequences of intergenic and intronic loci (hereafter, “non-coding loci”) were retrieved from those 32 plastomes. Loci smaller than 150 bp were excluded. We identified syntenic loci on the basis of their neighboring genes; i.e., loci that are flanked by the same genes/exons in the 32 sampled plastomes were identified as syntenic. Simple sequence repeat (SSR) elements located in the syntenic loci were detected using GMATo with the criteria that the “min length” for mono-nucleotide and multi-nucleotide SSRs was set to 8 and 5 units, respectively (Wang et al., 2013). Then we counted all of the SSR elements of syntenic loci for each genus, and the polymorphic SSR loci were counted one time.

To assess the sequence variability (SV) between congeneric species, the syntenic loci we identified were used for estimation and comparison. In total, we compared 84 pairs of congeneric species, including 1, 15, 55, 3, 1, 3, 6, 1, 1, and 1 pairs within Bletilla, Corallorhiza, Cymbidium, Dendrobium, Masdevallia, Phalaenopsis, Goodyera, Cypripedium, Paphiopedilum, and Vanilla, respectively. The sequences compared between each pair were aligned using MUSCLE 3.8.31 (Edgar, 2004) with the “refining” option. Pairwise nucleotide mutations and indel events were counted using DnaSP v5 (Librado and Rozas, 2009), with the exclusion of indels at the 5′- and 3′-ends of alignments. We adopted the methods of Shaw et al. (2005) and Downie and Jansen (2015) for calculating the SV of the examined loci. The formula was as follows: SV = (number of nucleotide mutations + the number of indel events)/(number of conserved sites + the number of nucleotide mutations + the number of indel events) × 100%.

In orchids, the loss of plastid ndh genes has independently occurred among different genera (Lin et al., 2015). The truncation or complete loss of ndh genes can lead to difficulties in alignment, caused by the high rates of nucleotide substitution or the presence of many deletions and insertions. Therefore, we only extracted sequences of 67 plastid protein-coding genes common to 53 orchids and eight other monocots (Supplementary Table S1) (without ndh genes). Sequence alignments were performed using MUSCLE. After removing all gaps and ambiguous sites, the alignments of the 67 genes were concatenated using SequenceMatrix 1.8 (Vaidya et al., 2011). The partitions of the 67 concatenated genes were analyzed using PartitionFinder v1.1.1 (Lanfear et al., 2012) and the result was incorporated into construction of maximum likelihood (ML) and Bayesian inference (BI) trees using RAxML 8.0.2 (Stamatakis, 2014) and MrBayes 3.2 (Ronquist et al., 2012), respectively. Fritillaria taipaiensis and Lilium longiflorum were designated as outgroups. The robustness of the ML tree was estimated using 1,000 bootstrap replicates. We performed two independent MCMC runs in analyses of the BI tree. Each run yielded one million generations. We collected one tree per 1,000 generations. The initial 25% of the collected trees were discarded as burn-in, and the remaining ones were used to estimate posterior probability.

The Illumina paired-end sequencing produced approximately 3.7 Gb of 73 bp pair-end reads for each species. The de novo assembly included 39,055 contigs for D. moniliforme, 44,793 contigs for G. schlechtendaliana, 10,613 contigs for P. armeniacum, and 15,676 contigs for V. aphylla. After comparisons were conducted using the plastomic sequences of G. fumata (NC_026773), 38 contigs were obtained with E values < 10^-10 and mean coverage depth > 30× for D. moniliforme, 83 contigs for G. schlechtendaliana, 87 contigs for P. armeniacum, and 47 contigs for V. aphylla. Two contigs (length = 86,931 bp and 35,926 bp) resulted in a nearly complete draft genome for D. moniliforme. Four contigs longer than 20 kb were used for G. schlechtendaliana assembly, and five contigs (longer than 20 kb) and three contigs (longer than 40 kb) were used to assemble the plastomes of P. armeniacum and V. aphylla, respectively. After assembly and gap closure, four complete plastomes were obtained.

The newly sequenced plastomes of D. moniliforme (Epidendroideae), G. schlechtendaliana (Orchidoideae), P. armeniacum (Cypripedioideae), and V. aphylla (Vanilloideae) are circular with IRs (i.e., IRA and IRB) separated by large single-copy and SSC regions (Figure 2). They range in size from 148,778 to 162,835 bp, with a GC content between 35.02 and 37.54% (Table 1). These four orchid plastomes differ in IR lengths, from 25,983 to 33,641 bp. The genes located in the IRs also vary greatly. The IRs of both D. moniliforme and G. schlechtendaliana contain a small fraction of the gene ycf1. By contrast, the IR of P. armeniacum encompasses three genes, ycf1, rps15, and psaC. That of V. aphylla includes ycf1, rps15, trnL, and a partial segment of ccsA, although the ycf1 has mutated with numerous stop codons (Figure 2).

Among the four orchid plastomes sequenced, only G. schlechtendaliana contains a full set of 11 plastid ndh genes. By contrast, D. moniliforme and P. armeniacum have truncated or lost 10 of them (only the ndhB gene is functional; Figure 2). No functional ndh gene was detected in the plastome of V. aphylla (Figure 2). Because the four sequenced orchids belong to four different subfamilies, their distinct loss/retention of ndh genes reflects multiple events of plastid gene loss during their evolution. Furthermore, the loss of plastid ndh genes and the expansion of IRs have caused drastic reductions in SSC regions. For example, the length of the SSC region in V. aphylla is only about one-eighth of that in G. schlechtendaliana (Figure 2).

Compared to the plastomes, which have full sets of ndh genes (Supplementary Table S2), the IR/SSC boundaries of Vanilla and Paphiopedilum have expanded to a highly anomalous position. We determined the degree of IR expansion/contraction based on the length of a region from the 5′ end of the ycf1 gene to the IR and SSC junction. A two-sided Mann–Whitney test showed that the IR lengths in Vanilla and Paphiopedilum plastomes were significantly different from plastomes with 11 functional ndh genes (P < 0.05). These results indicate that ndh genes play an important role in IR/SSC junction stability.

We identified 68 syntenic non-coding loci. SSRs within these loci were also counted (Supplementary Table S3). In these plastomes, the ratios of SSR-containing loci to SSR-lacking loci are significantly (two-sided Fisher’s exact test, all P < 0.05) higher in the SC than in the IR regions, indicating that the distribution of SSRs is dependent on their locations in plastomes. The mean values of pairwise interspecific SV were compared and the results are shown in Supplementary Table S3. The 68 syntenic intergenic spacers and introns were sorted into SC and IR loci, depending on their location. Our analysis, using a two-sided Mann–Whitney test, indicated that in all of the 10 orchid genera examined, the SC loci have significantly greater SV than the IR loci (all P < 0.05) (Supplementary Table S3). Moreover, mean SV was inversely correlated with mean GC content in all of the 10 orchid genera (Spearman’s r = -0.496 to -0.80, all P < 0.01), suggesting that highly mutated loci have also evolved toward the enrichment of AT nucleotides. To determine whether the evolution of SV is also conserved among the orchid genera, we conducted a correlation test of the 68 syntenic loci between and within different subfamilies. Except for Cymbidium vs. Bletilla (Spearman’s r = 0.767, P < 0.01), the SV values of the 68 syntenic loci within Epidendroideae were statistically uncorrelated (Corallorhiza vs. Phalaenopsis, r = 0.238, P > 0.05), slightly correlated (r = 0.238 to 0.488, all P < 0.05), or intermediately correlated (r = 0.513 to 0.675, all P < 0.05). Moreover, no strongly correlated relationship between different subfamilies was detected (Corallorhiza vs. Vanilla, r = 0.086, P < 0.05, others r = 0.243 to 0.675, all P < 0.05). These results suggest that in orchid plastomes, the evolution of SV is genus specific.

Figure 3 shows the top 10 loci that have the greatest SV for each genus. These loci are exclusively mutational hotspots. Obviously, most of them contain SSRs. However, none of the hotspots is common to all 10 genera, and only four of them (5′trnK-rps16, trnS-trnG, clpP-psbB, and rps16-trnQ) are present in more than four examined genera of Epidendroideae. In addition, among these four hotspots, only two, clpP-psbB and rps16-trnQ, are present in the two genera of Cypripedioideae. Therefore, we conclude that plastomic mutational hotspots are diverse among orchid genera.

We performed phylogenetic analyses of 67 concatenated plastid genes with or without partitions of sequences. Both ML and BI trees recovered a monophyly of the five orchid subfamilies, irrespective of whether or not the partitions of sequences were incorporated (Figure 4 and Supplementary Table S3). However, without partitioning, the ML and BI trees were inconsistent in the inferred positions of Neottia and Epipogium. As shown in Supplementary Figure S1, the ML tree placed both Neottia and Epipogium within Epidendroideae, whereas these two genera are placed as sister to Orchidoideae rather than to other species of Epidendroideae in the BI tree. However, with partitioning (Supplementary Table S4), both trees congruently placed Neottia and Epipogium within Epidendroideae (Figure 4 and Supplementary Table S4) and the trees no longer contradicted each other in their topologies (Figure 4 and Supplementary Table S1).

The Apostasioideae, whose members have multiple stamens, is the earliest diverging subfamily of orchids. The Vanilloideae, although consisting of one-stamened members, is closely related to a multi-stamened subfamily (i.e., Cypripedioideae) rather than to the other two one-stamened subfamilies (i.e., Epidendroideae and Orchidoideae). Therefore, on the basis of our phylogenetic analyses, the one-stamened orchids are paraphyletic. This suggests that the number of stamens evolved independently at least twice during the evolution of orchids.

In orchids, the lost plastid ndh genes were previously hypothesized to be functionally replaced by their putative nuclear homologs of plastid origin (Chang et al., 2006). However, analyses of nuclear transcriptome have failed to detect the complete set of the 11 ndh genes in some orchid species whose plastids lack ndh genes (Johnson et al., 2012; Su et al., 2013; Tsai et al., 2013; Lin et al., 2015). Considering the notable variation in the loss/retention of ndh genes among the four newly elucidated plastomes (Figure 2), we suggest that multiple ndh loss events have occurred during the evolution of orchids. This is in good agreement with the proposition that loss of ndh genes has occurred independently among orchid genera and subfamilies (Kim H.T. et al., 2015; Lin et al., 2015).

Plastid ndh proteins form a complex that joins cyclic electron transport at photosystem I (Martín and Sabater, 2010). Based on the existence of the PGR5-dependent cyclic electron transport pathway and the fact that no deleterious effects have been observed in ndh-deficient mutants growing under favorable conditions, Ruhlman et al. (2015) deduced that plastid ndh genes might be dispensable in contemporary plants. Nevertheless, the loss of ndh genes has led to structural changes in the plastomes of orchids, as shown by Kim H.T. et al. (2015) and us here.

The structure and gene content of the mycoheterotrophy orchid plastome are highly variable in photosynthetic orchids (Delannoy et al., 2011; Logacheva et al., 2011; Barrett and Davis, 2012). Our comparative analyses found that the Vanilla and Paphiopedilum plastomes have also expanded their IRs to a highly anomalous position. However, the mechanisms that underlie shifts of their IR boundaries are not clear. Kim H.T. et al. (2015) proposed that the instability of the IR/SSC junctions in orchids is strongly correlated with the deletion of the ndhF gene. Our two-sided Mann–Whitney test also indicated that ndh genes play an important role in IR/SSC junction stability. The loss of the full set of 11 plastid ndh genes has been reported in diverse photosynthetic seed plants, such as gnetophytes (Braukmann et al., 2009; Wu et al., 2009, 2011), orchids (Chang et al., 2006; Wu et al., 2010; Pan et al., 2012; Lin et al., 2015), pines (Braukmann et al., 2009; Wu et al., 2013), slender naiads (Peredo et al., 2013), and saguaros (Sanderson et al., 2015). With the exception of pines and saguaros, which have only one IR copy, all lineages mentioned above have expanded their IRs to include some genes (i.e., ycf1, rpl32, trnL, rps15, psaC, or ccsA) that are usually situated in the SSC regions of seed plant plastomes. Given that plastid ndh genes were helpful for the instability of IR/SSC junctions and that the expansion of IRs is common to ndh-deleted plastomes, we infer that shifts of IR boundaries might be associated with the variable loss of ndh genes in the Vanilla and Paphiopedilum plastomes. As the expansion of IRs directly causes gene duplication (Figure 2), whether it can benefit adaptation is worthy of further investigation.

It has been shown that plastome-wide comparisons facilitate the screening of mutational hotspots used for intraspecies discrimination (Ahmed et al., 2013; Hsu et al., 2014) and phylogenetic studies at the species level (Shaw et al., 2014; Downie and Jansen, 2015). In orchids, numerous non-coding plastid loci have been shown to be useful in species-level phylogenetic studies. For instance, the non-coding loci rpl32-trnL, trnE-trnT, trnH-psbA, trnK-rps16, and trnT-trnL may be markers for identifying species of Cymbidium (Yang et al., 2013), and trnS-trnG, psaC-ndhE, clpP-psbB, rpl16 intron, rpoB-trnC, trnT-psbD, rbcL-accD, rpl32-trnL, ccsA-ndhD, and ndhC-trnV may be useful for Phalaenopsis (Shaw et al., 2014). However, none of these was detected in IRs of the orchid plastomes, nor were the loci we report in Figure 3. These data reinforce the view that plastid substitution rates are considerably lower in IRs than in SC regions (Wolfe et al., 1987; Smith, 2015; Wu and Chaw, 2015).

Different mutational hotspots have been used for phylogenetic and identification analyses in orchid species, however, these have not been proposed to be diverse mutational hotspots shared among orchid genera (e.g., Neubig et al., 2009; Yang et al., 2013; Luo et al., 2014). In this study, our correlation analyses revealed that SV between and within different orchid subfamilies are uncorrelated or not strongly correlated, which suggests that in orchid plastomes, the evolution of SV is genus specific (Supplementary Table S3). Furthermore, the top 10 loci we screened as most likely containing the highest degrees of genetic variability in orchids are quite diversified (Figure 3). These results indicate that a weakly evolving locus in one orchid genus might be a quite variable locus in another genus.

Although highly variable mutational hotspots are diverse among orchid genera, we propose that the three loci 5′trnK-rps16, trnS-trnG, and rps16-trnQ might be powerful markers for genera within Epidendroideae, and clpP-psbB and rps16-trnQ might be useful for Cypripedioideae, for two reasons. The three hotspots 5′trnK-rps16, trnS-trnG, and rps16-trnQ are present in more than four genera of Epidendroideae, and the hotspots clpP-psbB and rps16-trnQ are present in the two examined genera of Cypripedioideae (Figure 3). Additionally, among these hotspots, 5′trnK-rps16, trnS-trnG, and rps16-trnQ are located in one of the three most variable plastome regions, matK to 3′trnG, as determined by Shaw et al. (2014) based on comparisons of non-coding loci among different plant genera. Our results are in good agreement with Shaw et al. (2014) and Downie and Jansen (2015). These authors identified the same highly variable loci in several disparate plant lineages.

Our findings will be helpful in identifying mutational hotspots representing the orchid subfamily level. However, there is still an indispensable need for the careful discovery and characterization of loci specific to all genera in orchids.

However, molecular phylogenetic studies of orchid to date have failed to agree on the placement of Cypripedioideae and Vanilloideae (Cameron et al., 1999; Freudenstein and Chase, 2001; Cameron, 2004, 2007; Górniak et al., 2010; Givnish et al., 2015; Kim H.T. et al., 2015). Recently, Givnish et al. (2015) and Kim H.T. et al. (2015) reconstructed ML trees of 39 and 14 orchid genera, respectively, using the concatenated nucleotide sequences of plastid genes. The results of our tree-based partitioned analyses (Figure 4) are congruent with theirs in the relationships among the five orchid subfamilies and in the monophyly of the clade Epidendroideae–Orchidoideae–Cypripedioideae. The monophyly of Epidendroideae, Orchidoideae, and Cypripedioideae, which is strongly supported in our plastid phylogenomics analyses, is also consistent with a number of previous studies that used increased taxon sampling or different molecular markers (e.g., Cozzolino and Widmer, 2005; Górniak et al., 2010; Lin et al., 2015). Neither their nor our results support the view suggested by morphological studies (Vermeulen, 1966; Rasmussen, 1985; Szlachetko, 1995) that the one-stamened subfamilies Epidendroideae, Orchidoideae, and Vanilloideae are monophyletic. Moreover, morphological studies have also shown that, in development, the single fertile anthers of the Vanilloideae are not homologous to those of the Orchidoideae and Epidendroideae (Freudenstein et al., 2002). Therefore, we conclude that one-stamened orchids are paraphyletic.

The generic relationships found in our partitioned analyses within Epidendroideae are largely congruent with those of recent studies (e.g., Górniak et al., 2010; Freudenstein and Chase, 2015; Givnish et al., 2015; Kim H.T. et al., 2015), with only some genera being weakly supported. This may be attributable to the absence of numerous photosynthesis-related genes from Corallorhiza, Epipogium, and Neottia, resulting in large amounts of missing data for these taxa (Kim H.T. et al., 2015). Moreover, incomplete taxon sampling has likely also resulted in the inconsistent placement of the tribe of Epidendreae and Sobralieae found in previous studies (Górniak et al., 2010; Freudenstein and Chase, 2015; Givnish et al., 2015) and ours. We anticipate that increases in the breadth of taxon sampling will improve the resolution of orchid phylogenies.