Phylogeny and Historical Biogeography of Paphiopedilum Pfitzer (Orchidaceae) Based on Nuclear and Plastid DNA

The phylogeny and biogeography of the genus Paphiopedilum were evaluated by using phylogenetic trees derived from analysis of nuclear ribosomal internal transcribed spacer (ITS) sequences, the plastid trnL intron, the trnL-F spacer, and the atpB-rbcL spacer. This genus was divided into three subgenera: Parvisepalum, Brachypetalum, and Paphiopedilum. Each of them is monophyletic with high bootstrap supports according to the highly resolved phylogenetic tree reconstructed by combined sequences. There are five sections within the subgenus Paphiopedilum, including Coryopedilum, Pardalopetalum, Cochlopetalum, Paphiopedilum, and Barbata. The subgenus Parvisepalum is phylogenetic basal, which suggesting that Parvisepalum is comprising more ancestral characters than other subgenera. The evolutionary trend of genus Paphiopedilum was deduced based on the maximum likelihood (ML) tree and Bayesian Evolutionary Analysis Sampling Trees (BEAST). Reconstruct Ancestral State in Phylogenies (RASP) analyses based on the combined sequence data. The biogeographic analysis indicates that Paphiopedilum species were firstly derived in Southern China and Southeast Asia, subsequently dispersed into the Southeast Asian archipelagoes. The subgenera Paphiopedilum was likely derived after these historical dispersals and vicariance events. Our research reveals the relevance of the differentiation of Paphiopedilum in Southeast Asia and geological history. Moreover, the biogeographic analysis explains that the significant evolutionary hotspots of these orchids in the Sundaland and Wallacea might be attributed to repeated migration and isolation events between the south-eastern Asia mainland and the Sunda Super Islands.


INTRODUCTION
Paphiopedilum is a genus of tropical Asiatic origin, and its range extends eastward, reaching the Philippines, Southeast Asia, Borneo, and the Malay Archipelago, crossing Wallace's Line into Sulawesi, the Moluccas, New Guinea, and the Solomon Islands (Cribb, 1998). Tracking back to the geological history of Southeast Asian, the Palawan, Mindoro, Zamboanga, and the adjacent small islands are the older islands of the Southern Philippines. These regions are located on the border of the Eurasian Plate and have been shifting away from the mainland mass by tectonic collision since the early Miocene (~30 Mya) and the shell of the older plate was merged to Borneo until 5~10 Mya (Karig et al., 1986;Stephan et al., 1986;Hall, 1996). In contrast, most of the Philippine islands formed less than 5 Mya (Aurelio et al., 1991;Quebral et al., 1994). In addition, the Sundaland was comprised of the Malay Peninsula, Sumatra, Java, Sulawesi, and Borneo and merged with Bali, the Philippines, and even New Guinea/Australia into Sunda Superland interconnecting by land bridge during the last glacial period (0.01~1.8 Mya) (van Oosterzee, 1997). Since the last glacial period, species migrated forward and backwards between these regions and isolated after the last glacial maximum (LGM), causing the broken of Sunda Superland (Tsai et al., 2015).
The chloroplast primers for the atpB-rbcL, trnL-trnF spacer, and trnL intron are useful for phylogenetic studies at the intrageneric level. The primers for the trnL-trnF spacer and trnL intron developed by Taberlet et al. (1991) have been applied for inferring phylogenies at the intrageneric level (Goldblatt et al., 2002;Hodkinson et al., 2002;Mogensen, 1996;Van Raamsdonk et al., 2003), and have also been used successfully on Orchidaceae (Tsai et al., 2012). The atpB-rbcL regions are high length differences due to frequent occurrence of indels and are often used in combination with other primers to provide more information (Yoshinaga et al., 1992;Chiang et al., 1998;Von Konrat et al., 2010). Therefore, this study aims to further elucidate the phylogeny of Paphiopedilum through analysis ITS (internal transcribed spacer) sequences and three non-coding plastid DNA sequences (trnL intron, trnL-F, and atpB-rbcL spacers). In addition, the biogeography of this genus is clarified based on the phylogenetic tree derived from the molecular evidence.

Plant Materials
Seventy-eight taxa of Paphiopedilum and two outgroups from genus Phragmipedium were used in this study ( Table 1). All leaf materials were taken from living plants in the greenhouse of the Kaohsiung District Agricultural Improvement Station (KDAIS) in Taiwan.

DNA Extraction, PCR Amplification, and Sequencing
Total DNA was extracted from fresh etiolated leaves by using the cetyltrimethylammonium bromide (CTAB) method (Doyle and Doyle, 1987). Approximate DNA yields were determined by using the spectrophotometer (model U-2001, Hitachi).
The PCR reaction was used to amplify nuclear ribosomal ITS sequence and chloroplast (cp) DNA fragments trnL intron and the trnL-trnF spacer, atpB-rbcL spacer. ITS primers were designed from conserved regions of the 3' end of the 18S rRNA gene and the 5'end of the 26S rRNA gene using sequences from different species present in GenBank. Universal primers for trnL intron and the trnL-trnF spacer were referenced from Taberlet et al. (1991). Primer sequences for amplifying of the atpB-rbcL spacer were designed from the conserved regions of the 3' end of the atpB gene and the 5'end of the rbcL gene of chloroplast DNA using sequences of different species obtained from GenBank. Detailed amplification conditions and primer sequences are given in Supplementary Table S1. All PCR products were separated by agarose gel electrophoresis (1.0%, w/v in TBE) and were recovered using glassmilk (BIO 101,California).
PCR products were directly sequenced using the dideoxy chain-termination method on an ABI377 automated sequencer with the Ready Reaction Kit (PE Biosystems, California) of the BigDye™ Terminator Cycle Sequencing. The PCR reaction primer sequences were used as sequencing primers. Each sample was sequenced two or three times to confirm the sequences. Reactions were performed as recommended by the product manufacturers.

Sequence Alignment and Phylogenetic Reconstruction
The sequence alignment was determined using the ClustalW multiple alignment program in BioEdit (Hall, 1999), and four regions were combined for the following analysis. The alignment was checked, and apparent alignment errors were corrected by hand. Indels (insertion/deletions) were treated as missing data. For phylogenetic reconstruction, two Phragmipedium taxa treating as outgroups were sequenced to resolve whether all ingroup taxa formed a monophyletic lineage. The best-fitting substitution model was selected (Supplementary Table S2) by a model test using MEGA 6.0 (Tamura et al., 2013). Tamura 3parameter model (T92) using a discrete Gamma distribution (+G) was selected for following neighbor-joining (NJ) phylogenetic reconstruction. The general time reversible (GTR) using a discrete gamma distribution (+G) and considering the proportion of invariable sites (+I) were chosen for following divergence time estimation using the Yule model methods in BEAST 1.8.0 (Drummond and Rambaut, 2007;Drummond et al., 2012). The phylogenetic tree for the combined multiple sequence datasets used equally weighted characters. Moreover, because the sequence data of the four genera (Mexipedium, Selenipedium, Cypripedium, and Goodyera) in NCBI is limited, only two sets of fragment data (ITS: Mexipedium xerophyticum-MK161260.1; Selenipedium aequinoctiale-JF825977.1; Cypripedium_macranthos-KT338684.1; Goodyera_procera-MK451741.1and trnL-trnF spacer: Mexipedium xerophyticum-FR851215.1; Selenipedium aequinoctiale-JF825975.1; Cypripedium_macranthos-JF797026.1; Goodyera_procera-MK451782.1) are used as an additional analysis and compared with the data using only genus Phragmipedium as outgroup. The results of six outgroups are showed in Supplementary Data. 1 | Names of specimens, geographical distribution, source, and GenBank accession numbers for sequences of the internal transcribed spacer (ITS) of ribosomal DNA (rDNA), the plastid trnL intron, the trnL-F spacer, and the atpB-rbcL spacer.  Genetic relationships were determined using NJ in the MEGA 6.0 (Tamura et al., 2013), maximum parsimony (MP) in PHYLIP 3.68 (Felsenstein, 2004), and maximum likelihood (ML) in MEGA 6.0 (Tamura et al., 2013). Bootstrapping (1,000 replicates) was carried out to estimate the support for NJ, MP, and ML topologies (Felsenstein, 1985;Hillis and Bull, 1993). The strict consensus parsimonious tree was then constructed by using the MEGA 6.0 (Tamura et al., 2013).

Divergence Time Estimation
The combined chloroplast DNA (cpDNA) dataset was used to estimate the divergence times using the Bayesian Yule model methods (BEAST version 1.7.5). The characteristic of uniparental inheritance in cpDNA prevents the interference of recombination introgression on phylogenetic reconstruction (Drummond et al., 2012). The general-time reversible (GTR) model with estimates of invariant sites (+I) and gamma-distributed among site rate variation (+G) in all matrices without partitions model was determined by the nucleotide substitution model test, conducted in MEGA 6.0 (Tamura et al., 2013).
To estimate the divergence time, two strategies, the strict and relaxed clock models, were adopted. For the relaxed clock, the calibration point at the most recent common ancestor (MRCA) of Paphiopedilum and Phragmipedium for 22 Ma (Gustafsson et al., 2010) were used to calculate the divergence times of each node. However, since there is no suitable fossil record to correct the calibration of divergence times for the ingroup, we recalculate the divergence time by strict clock model for consistency. For the strict clock, the mean substitution rate of 1.82 × 10 −9 subs/ site/year with the lower and upper limits 1.11 × 10 −9 subs/site/year and 2.53 × 10 −9 subs/site/year, respectively, were used for the cpDNA spacers in Phalaenopsis amabilis complex (Tsai et al., 2015).
We conducted four independent runs of a Yule prior and four Markov Chain Monte Carlo (MCMC) chains with a different starting seed. The first 10% simulations were discarded (burn-in) in a total of 10 8 generations. For thinning, one tree was reserved every 10,000 trees, and finally, 10,000 trees were left to calculate the posterior probability of each node. The effective sample size (ESS) > 200 was used as a criterion to check whether the sampling (simulations) is proper and is reaching a stationary distribution by Tracer v1.6 (Rambaut et al., 2018). Four independent-runs results, including the log file and tree file, were combined with the assistance of LogCombiner 1.6.1 (Drummond and Rambaut, 2007). TreeAnnotator 1.6.1 (Drummond and Rambaut, 2007) was used to summarize a consensus tree with a criterion of the maximum clade reliability using the mean heights option. The final consensus tree was drawn by FigTree 1.3.1 (Rambaut, 2009).

Biogeographic Inference Using Reconstruct Ancestral State in Phylogenies
The Statistical dispersal-vicariance analysis was used to assess the biogeographic patterns of Paphiopedilum species [Statistical Dispersal-Vicariance Analysis (S-DIVA), (Yu et al., 2010)]. Bayes-Lagrange Statistical dispersal-extinction-cladogenesis (S-DEC) model (Ree and Smith, 2008) was performed in Reconstruct Ancestral State in Phylogenies (RASP) 3.2 (Yu et al., 2015) to distinguish the events of vicariance, dispersal, and extinction. Five geographic areas were determined mainly according to Myers et al. (2000) with a little modification to illuminate the vicariance, dispersal and extinction events of Paphiopedilum species. The hotspot areas in South-Central China and Indo-Burma with the Malay Peninsula were combined as area A, consisting of China, Nepal, India, Bhutan, Burma, Thailand, Malaysia, Cambodia, Vietnam, and Laos. The hotspot "Sundaland" including Borneo, Java, and Sumatra were set as the area B. We move the Malay Peninsula from the area "Sundaland" to area A due to the integrity of the current landmass. The hotspots "Wallacea" (include Sulawesi and Moluccas) and "Philippines" were set as area C and area E. Respectively, islands of New Guinea eastern from the Wallacea are defined as area D. Species of outgroup were all defined as the area I. These two outgroup species are distributed in Ecuador, Peru, Costa Rica, Panama, and Colombia. The real distribution of outgroups is too far from the areas of the species in this study. Therefore, the ranges of outgroups are assigned to a new area in which none of the ingroup species occurs (Yu et al., 2014). The ML tree topologies were used in S-DIVA analysis.

Sequence Alignment and Characteristics
The lengths of the ITS sequences obtained from the Paphiopedilum and outgroup samples were similar to those reported for a broad example of angiosperms (Baldwin, 1992;Baldwin et al., 1995). The alignment length of the ITS sequence is 735 nucleotides, of which 343 were identified as variable sites with 235 potentially parsimony informative sites. The average genetic distance between the 78 Paphiopedilum samples was 0.039 in ITS, and the average genetic distance between the 78 Paphiopedilum species was 0.01 in cpDNA. The alignment of combined plastid DNA fragments contained a total of 2,409 characters, of which 872 were identified as variable sites with 588 potentially parsimony informative sites. Since the sequences of three samples of every species are the same within species, only one sequence per species was used for the analyses and deposited in NCBI GenBank. The accession numbers of the nuclear ribosomal ITS sequences and the three fragments of plastid DNA from the 78 Paphiopedilum taxa and the two outgroup samples (from the genus Phragmipedium) are shown in Table 1.

Phylogeny Reconstruction
Both NJ and MP trees revealed a monophyletic relationship of 78 Paphiopedilum taxa with high bootstrap supports ( Figure 2). Moreover, the use of six outgroups for phylogeny reconstruction showed similar bootstrap supports (Supplementary Figure S1).  Figure S1)] for the monophyly.
In addition, the molecular data demonstrates that a newly described variety, P. micranthum var. eburneum, is closely related to P. malipoense based on the plastid DNA within subgenus Parvisepalum, which is inconsistent with the inference by nuclear ITS and combined data. In ITS tree, P. micranthum var. eburneum is sister with P. micranthum var. micranthum (Figures 2 and 3). Therefore, we infer a hybridization event between the ancestor of P. micranthum and P. malipoense that lead to a plastid capture in P. micranthum var. eburneum.

Demographic History and Historical Biogeography Inference
Complicated evolutionary processes of continuous and episodic dispersal, vicariance, and extinctions determined the current geographic distribution of genus Paphiopedilum. Since the most probable ancestral areas located on continental Asia (area A in Figure 5), dispersal events seem to determine the extant distributions of subgenera largely. The results supported vicariance events on nodes 159, 146, and 109 shown in Table 2 and Figure 5, and on nodes 163, and 132 (Supplementary Table S3 and Figure S3). The node 147 revealed dispersal events among section Coryopedilum/Pardalopetalum of Paphiopedilum and other sections of Paphiopedilum causing by migration route from area A (China, Nepal, India, Bhutan, Burma, Thailand, Malaysia, Cambodia, Vietnam, and Laos) to B area (Sumatra, Borneo, and Java). Meanwhile, the nodes 145, 129, 114, and 102 also revealed dispersal events from north to south, according to Sundaland and Sunda Super Islands.
Furthermore, the use of six outgroups for dynamic historical inference showed similar patterns (Supplementary Table S3 and Figure S3). Only the nodes 146 and 109 were detected vicariance event causing by the geological separation between Indochina and Sumatra/Borneo/Java. In addition, in subgenus Paphiopedilum, 2 vicariance and 10 dispersal events were detected, which suggesting a significant dispersal process affected biogeographical patterns in shaping the current distribution in the subgenus Paphiopedilum. Areas A and B might be the two possible ancestral areas and likely shaped by several complicated dispersal events in the subgenus Paphiopedilum.

Systematics Revision of Genus Paphiopedilum
In general, our phylogenetic inference is mostly congruent with that of Cox et al. (1997), Cribb (1997b), and Guo et al. (2015). In the genus Paphiopedilum, tessellated leaves, single flowers with broad elliptic to subcircular petals, and a sizeable thin-textured lip characterize subgenus Parvisepalum in southwest China and Vietnam (Cribb, 1998). Within this subgenus, the phylogenetic topography and the divergence time of at least 4.30 Mya rejected the previous hypothesis of the sister-species relationship between P. armeniacum and P. delenatii (Cribb, 1983) (Figure 4). The geographical distribution of these two species is also separated (Yunnan, China for P. armeniacum, and Vietnam for P. delenatii) (Cribb, 1998). Furthermore, a newly described variety, P. micranthum var. eburneum, is phylogenetically close to P. malipoense in maternalinherited plastid DNA but close to P. micranthum in biparentalinherited nuclear ITS sequences ( Figure 5), suggesting that P. micranthum var. eburneum is a natural hybrid between the maternal parent P. malipoense and the paternal parent P. micranthum and experienced the event of chloroplast FIGURE 2 | Phylogenetic relationships using neighbor-joining (NJ), maximum parsimony (MP), and maximum likelihood (ML) resulting from analysis of the combined data matrix (nuclear ribosomal ITS, plastid trnL intron, trnL-F spacer, and atpB-rbcL spacer) from 78 Paphiopedilum and 2 outgroup species. Only the strict consensus of all most parsimonious trees (MP trees) are showed in this figure, and the bootstrap values > 50% are shown on each branch for NJ/MP/ML between major lineage.
capture. The overlap of the geographical distribution of these three taxa also supports this hypothesis (Cribb, 1998). In addition, ITS sequences are usually concertedly evolved via unequal crossing-over (Schlotterer and Tautz, 1994) and biased gene conversion (Hillis et al., 1991), which results in sequence homogeneity between paralogs (Maynard, 1989).
The monophyly of subgenus Brachypetalum inferred in this study is congruent with the inference by Cox et al. (1997). The subgenus Brachypetalum is geographically confined to Southeast Asia (Cribb, 1998). Albeit overlapping distribution with subgenus Parvisepalum (Cribb, 1998), subgenus Brachypetalum is phylogenetically separated, consistent with the distinguishable leaf anatomy between these two subgenera (Cribb, 1998). Both molecular and morphological evidences support the independent taxonomic treatment between subgenera Brachypetalum and Parvisepalum (Karasawa, 1982;Karasawa and Saito, 1982;Cribb, 1997b), but object with Atwood's (1984) opinion of taking the subgenus Parvisepalum as a synonym of Brachypetalum.
FIGURE 4 | Results of calescence time estimations performed with BEAST 1.8.0 for the from 78 Paphiopedilum taxa based on the combined data matrix (nuclear ribosomal ITS, plastid trnL intron, trnL-F spacer, and atpB-rbcL spacer). The black numbers in each node are using the mean rate of 1.82 × 10−9 subs/site/year, with the lower and upper limits 1.11 × 10−9 subs/site/year and 2.53 × 10−9 subs/site/year (Tsai et al., 2015). The red numbers in each node are using the Gustafsson et al. (2010) fossil calibration time data to calibrate the divergence time. The monophyletic subgenus Paphiopedilum can be morphologically and phylogenetically subdivided into five sections: Coryopedilum, Pardalopetalum, Cochlopetalum, Paphiopedilum, and Barbata (Cox et al., 1997;Cribb, 1997b). Section Coryopedilum is characterized by its plain green, straplike leaves, and multi-flowered inflorescences, which flowers are blooming simultaneously (Cribb, 1998). This section distributes throughout Borneo, Sulawesi, New Guinea, and the Philippines (Cribb, 1998). Except placing Paphiopedilum parishii and Paphiopedilum dianthum into section Pardalopetalum from section Coryopedilum, Cribb (1997b) agreed with Atwood (1984) that section Pardalopetalum is independent from section Coryopedilum taxonomically, according to the ITS analysis (Cox et al., 1997) and similar green strap-like leaves and staminodes (Cribb, 1997b), which is consistent with our phylogenetic inference. However, the only character that separates sections of Pardalopetalum and Coryopedilum is the morphology of staminode. Whether this single character is sufficient to characterize them as separating sections should be re-evaluated with more evidence.
Unlike the simultaneous bloom of section Coryopedilum, section Cochlopetalum flower in succession, and their flowers bear elliptic bracts, linear, spirally twisted, spreading, ciliate petals, and a pot-shaped spotted lip (Cribb, 1998). Section  (Yu et al., 2010)] and Bayes-Lagrange Statistical dispersalextinction-cladogenesis (S-DEC) model (Ree and Smith, 2008) performed in Reconstruct Ancestral State in Phylogenies (RASP) 3.2 (Yu et al., 2015). Phylogenetic relationships of the 78 Paphiopedilum species, plus the two outgroups Phragmipedium pearcei and Ph. longifolium, obtained from sequence judgments of the Cochlopetalum distributes in Sumatra and Java only (Cribb, 1998). The extensive section Barbata is the sister of section Cochlopetalum, also characterized by a solitary flower with a lip and prominent incurved side-lobes, but the leaf tessellated (Cribb, 1998). The morphological dissimilarity and reciprocally monophyletic relationship indicate that, despite recently diverged, these two sections should be independent taxonomically.

Biogeography and Evolutionary Trends
The clade of genus Paphiopedilum is coalesced to 7.09 or 5.72 Mya, similar to the estimate of 7.62 Mya by Guo et al. (2015). The flower morphology of subgenus Parvisepalum is intermediate between other subgenera of Paphiopedilum and Cypripedium (Chen and Tsi, 1984), which could be explained by the earlier divergence of subgenus Parvisepalum in genus Paphiopedilum (Guo et al., 2012). Presently, the genus Cypripedium is distributed throughout worldwide temperate zones (Cox et al., 1997), with China as a center for species diversity (Cribb, 1997a). Therefore, genera Paphiopedilum and Cypripedium have most likely diverged in mainland Asia (Chen and Tsi, 1984).
However, genus Paphiopedilum was suggested as the sister with two American genera Phragmipedium and Mexipedium according to morphology, plastid rbcL (Albert, 1994), ITS (Cox et al., 1997), and both nuclear and plastid genes (Guo et al., 2012). These inferences are conflict to the hypothesis of the divergence between Paphiopedilum and Cypripedium in China, but implied the divergence of Paphiopedilum from the group of Phragmipedium + Mexipedium, by which the slipper orchids (Cypripedioideae) were hypothesized widespread throughout North America and Asia in the past (Atwood, 1984;Albert, 1994;Cox et al., 1997).
Subgenus Parvisepalum in southwest China and Vietnam diverged earlier from the other subgenera of genus Paphiopedilum. The coalescence time of subgenera Parvisepalum, Brachypetalum, and Paphiopedilum were tracking back to the Upper Miocene (Guo et al., 2015). Subgenus Brachypetalum in mainland Southeast Asia was descended from the subgenus Parvisepalum inferred by S-DIVA ( Figure 5), which agrees with other disjunctions at the Southern China and Indochina (Guo et al., 2015) or Sunda Shelf and New Guinea/Australia (Tougard, 2001;Lohman et al., 2011;Tsai et al., 2015).
The subgenus Paphiopedilum is further descended and evolved quickly in the Sunda Shelf. A land bridge might connect Mindoro, Palawan, Borneo, the Malay Peninsula, Borneo, Sumatra, Java, Bali, and various parts of the Philippines when sea levels falling during Pleistocene (about 0.01~1.8 Mya) (van Oosterzee, 1997) ( Figure 1). Surfaced land bridge connected these regions and was beneficial to the interisland and continent-island dispersal in Southeast and South Asia (van Oosterzee, 1997). The clade of Coryopedilum + Pardalopetalum was the first derived in subgenus Paphiopedilum based on the phylogenetic tree, which reflects in the sympatric distribution of subgenus Brachypetalum and section Pardalopetalum ( Figure 6). Following this clade formation, sections Paphiopedilum and Barbata were subdivided and dispersed throughout Southeast Asian archipelagoes across the land bridge during the glacials. The southward expansion from continental Asia into the Greater Sunda Islands through the Indochina and Malay Peninsulas were also reported in other taxa, e.g., Lithocarpus (Fagaceae) (Yang et al., 2018). Such colonization events between continental Asia and the Greater Sunda Islands mostly occurred during Miocene and Plio-Pleistocene (de Bruyn et al., 2014). As a "corridor," Indochina reveals high flora diversity and the high species richness, which facilitates the in situ speciation (de Bruyn et al., 2014).
Another flora diversity hotspot is Borneo (de Bruyn et al., 2014), which is also important for the genus Paphiopedilum. The section Cochlopetalum, which is found only in Sumatra and Java, might represent a group derived from Borneo. The S-DIVA inferred multiple times dispersal events sourced from Borneo with two vicariances to illustrate the current distribution of the five sections within subgenus Paphiopedilum ( Figure 5 and Table 2). The Borneo is the second original center of Paphiopedilum. The tropical forests and rugged topography harbor diversified niches opened to the speciation of organisms. The repeated submergence and emergence of land bridges could promote repeated genetic isolation and gene flow between closed related taxa. During the Plio-Pleistocene glacial oscillations. This process accelerates the diversification rates in the Sunda Super Islands.
Dispersal and vicariance events that exposed geographic isolation among taxa might be due to the land bridge submergence (Chiang et al., 2009;Chiang et al., 2013;Ge et al., 2012;Ge et al., 2015;Hsu et al., 2013) in Sunda Shelf and Sunda Super Island during the Pleistocene ( Figure 5 and Table 2) (Guo et al., 2015;Tsai et al., 2015). Cenozoic collision accompanied by a cyclical climate (glacial oscillations) caused by the fragmentation of the Sunda Super Islands (de Bruyn et al., 2014). The Borneo, Java, Sumatra, and the southern Philippines belong to the Sunda plate, Sulawesi is composed of broken plates, and the Moluccas and New Guinea belong to the Australian plate (Hall, 1998). The deep trenches between these plates cause segregation of species between Sundaland and the islands in the east. Because of the seed germination relies on symbiotic fungi, the geographical isolation maybe not only influences orchid itself but also in symbiotic fungi. However, these inferences still need further verification.

CONCLUSIONS
In summary, the origin and coalescence time of genus Paphiopedilum tracked back to Southern China/Eastern Indochina since late Miocene and early Pliocene, while the range expansion and species divergence were related to sea-level fluctuations during the Plio-Pleistocene glacial cycles. Historical geological barriers shaped a pattern of vicariance among disjunct distributed subgenera after isolated ancestral populations. The ancestral taxa of subgenus Paphiopedilum migrated from Southern China/Eastern Indochina to south which developed quickly in the Sunda Shelf. Due to the submergence of the Sunda Shelf and Sunda Super Island, species of subgenus Paphiopedilum dispersed with isolation between islands as well as subsequent in situ speciation within islands from other taxa within section or subgenus, which accelerated species divergence in subgenus Paphiopedilum. Paphiopedilum distributes in four of 25 biodiversity hotspots (Myers et al., 2000), the Indo-Burma, Sundaland, Wallacea, and Philippines, where are also the "major evolutionary hotspots" (de Bruyn et al., 2014). It is suggests that rich and fascinating historical biogeographic events have created rich species diversity there, such as the case of Paphiopedilum. However, deforestation has caused the so-called "empty forest syndrome" (de Bruyn et al., 2014). We hope that these areas will not become extinction hotspots, even though they are almost now.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
Conceived and designed the experiments: C-CT and Y-CC. Performed the experiments: C-CT, P-CL, Y-ZK, C-HC, and Y-CC. Analyzed the data: C-CT, P-CL, Y-ZK, C-HC, and Y-CC. Contributed reagents/materials/analysis tools: C-CT, P-CL, Y-ZK, C-HC, and Y-CC. Wrote the paper: C-CT, P-CL, and Y-CC. Conceived of the study, edited the manuscript, and approved the final manuscript: C-CT, P-CL, Y-ZK, C-HC, and Y-CC.

ACKNOWLEDGMENTS
The achievement of this study is dedicated to the memory of the scientific contributions of the first author C-CT, who died of a

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2020.00126/ full#supplementary-material SUPPLEMENTARY FIGURE S1 | Phylogenetic relationships using Maximum Likelihood resulting from analysis of the combined data matrix (nuclear ribosomal ITS, and trnL-F spacer) from 78 Paphiopedilum and 6 outgroup species.
SUPPLEMENTARY FIGURE S2 | Results of calescence time estimations performed with BEAST 1.8.0 for the from 78 Paphiopedilum taxa and 6 outgroup species based on the combined data matrix (nuclear ribosomal ITS, and trnL-F spacer).
SUPPLEMENTARY FIGURE S3 | Ancestral distributions reconstructed by the Statistical dispersal-vicariance analysis (S-DIVA, Yu et al., 2010) and Bayes-Lagrange Statistical dispersal-extinction-cladogenesis (S-DEC) model (Ree and Smith, 2008) performed in RASP 3.2 (Yu et al., 2015). Phylogenetic relationships of the 78 Paphiopedilum species, plus the six outgroups, obtained from sequence judgments of the combined sequence and generated by BEAST. The distribution areas of extant the 78 Paphiopedilum species species are marked in capitals A-E and I ((A) China, Nepal, India, Bhutan, Burma, Thailand, Malaysia, Cambodia, Vietnam, and Laos; (B) Sumatra, Borneo and Java; (C) Sulawesi and Moluccas; (D) New Guinea; (E) Philippines; and (I) outgroup), respectively. The green and blue circles indicate the vicariance and dispersal events obtained from the RASP analysis, respectively.