Previous phylogenetic analyses based on a small number of DNA loci or limited taxon sampling failed to robustly reconstruct the backbone of the Paris tree [21–24, 44–45]. Including all currently recognized species, the plastome analysis fully resolved interspecific relationships of Paris with strong support at most nodes. Our study further confirms that phylogenetic analysis based on more DNA loci with greatly increased number of phylogenetically informative characters can significantly improve resolution at low taxonomic levels [36, 43, 52, 53, 60].
The plastid analyses strongly support the monophyly of Paris and recovered five strongly supported major clades that correspond to the previously proposed sections by Ji et al. . Within Paris, successive divergence along the spine of the tree of sections Paris, Kinugasa, Axiparis, Thibeticae, Axiparis, and Euthyra was inferred. This divergence pattern can be supported by some morphological characters (Fig. 5). Briefly, the slender rhizome and round more berry-like fruit distinguish section Paris from the rest of the sections. Nevertheless, seeds without an enveloping sarcotesta (or aril, presumably a plesiomorphic character) separate sections Paris and Kinugasa from the rest. Although species of the thick rhizome clade (including sections Kinugasa, Thibeticae, Axiparis, and Euthyra) commonly have angular fruits, section Kinugasa is distinctive among these in possessing showy white sepals. The third diverging section Thibetica is similar to Euthyra in having dehiscent capsules, but its seed morphology (with an incomplete aril) is analogous to that of Axiparis. Therefore, the plastid tree is highly congruent with the morphologically based classification of Ji et al. .
The plastid tree provides valuable insights for resolving the long-standing disagreements in classification of Paris. All Paris species share the morphological synapomorphies of flowers and leaves being 4–15-merous compared with the trimerous condition of Trillium (Fig. 5), and monophyly of both genera was fully supported, making it reasonable to recognize Paris as a single genus [7, 17, 21] rather than dividing it into three genera (Daiswa, Kinugasa, and Paris s. s.) . Also, the new tree further supports the taxonomic treatment of Ji et al.  by combining sections Dunnianae, Fargesianae and Marmoratae with Euthyra. Given its economic importance, resolution of the long-standing taxonomic taxonomic disputes will be conducive to exploration and protection of Paris species.
The current taxonomy of P. polyphylla with five varieties is not supported by either the plastid or nrDNA results. The varieties should be probably recognized as distinct species, but it is also likely, given the cytonuclear discordance observed for these accessions (Fig. 2) that hybridization may be involved in their origins. Further study of this group with more appropriate population genetic and cytological techniques is warranted.
Similar to the previous study of Ji et al. , we detected discordance between the nrDNA and plastid trees for Paris at both deep and shallow nodes. Cytonuclear incongruence is a fairly common phenomenon in plant phylogenetics [25, 61–63]. In most cases, the nuclear tree is more congruent with morphological taxonomy [43, 56, 61, 62, 64, 65], and such incongruence can be mainly attributed to incomplete sorting of cytoplasmic polymorphisms or introgression of the cytoplasmic genome from one species into the nuclear background of another by hybridization [25, 63, 66, 67]. However, in our study the plastid tree of Paris is largely consistent with morphological evidence, suggesting nrDNA introgression without cytoplasmic gene flow could be responsible for the discordance detected [66–69].
Given that the discordance observed in Paris was likely due to phenomena affecting the nrDNA tree, which involved rapid gene conversion of one parental copy [70–75], the phylogenetic relationships recovered by this dataset may not be representative of that obtained with other parts of the nuclear genome not subject to gene conversion. Nonetheless, these results may provide useful information about past hybridization events that otherwise may not be the majority pattern in the nuclear genome. Failure to recover Paris as a monophyletic in the nrDNA tree (Fig. 2) suggests there could have been hybridization between section Paris and Trillium after section Paris split form the rest of the genus. It is noteworthy that the largest eukaryotic genome, that P. japonica [76, 77], was supposed to be an allo-octaploid between Paris and Trillium according to previous cytological investigations . This hybrid hypothesis is supported by its position as sister to section Paris and Trillium in the nrDNA tree. Likewise, the non-monophyly of sections Axiparis and Euthyra observed in the nrDNA tree (Fig. 2) can also be attributed to ancient hybridization. The sister relationship of P. luquanensis/P. mairei (section Euthyra) and the clade of four species of section Axiparis (P. dulongensis, P. forrestii, P. rugosa and P. tengchongensis) suggest that hybridization could have occurred between the ancestors of these taxa. Additionally, extensive discordance among species of section Euthyra (Fig. 2) supports the conclusions of the previous study that natural hybridization between species of section Euthyra is likely if the pollinators are the same, but little is known about this aspect of the biology of Paris. Experimental manual outcrossing has been effective between most of these species . Interspecific hybridization is the likely cause of the cytonuclear discordance observed between species in the section. Additionally, as mentioned above, there is discordance for the positions of the varieties of P. polyphylla, suggesting that hybridization may also have played a role in their origins.
Biogeography and lineage diversification
Because strong cytonuclear discordance was detected in Paris and the plastid tree agrees well with morphologically based classification, we address the biogeography and historical diversification of Paris based on the plastid dataset. The S-DIVA analysis recovered northeastern Asia and northern China as the ancestral area of Paris. Associated with a dispersal and a vicariance event (Fig. 1), the crown node of Paris was dated at 37.9 Mya (Fig. 3), near the Eocene/Oligocene boundary, when the global climatic deterioration  led to the expansion of vegetation adapted to drier and colder climates in large parts of Eurasia . Therefore, early divergence of Paris may have been driven by these events. Also, S-DIVA analysis revealed that the divergence of the Japanese endemic species, P. japonica and P. tetraphylla was triggered by two independent vicariance events (Fig. 1). Their divergence times of 27.63 (P. japonica) and 19.15 (P. tetraphylla) Mya, around the Oligocene/Miocene boundary corresponds to the opening of the Japan Sea, which separated the Japan islands from the continental East Asia .
In the thick rhizome clade, the S-DIVA analysis inferred three dispersal events (Fig. 1), which were dated at 18.64, 13.17 and 8.23 Mya (Fig. 3). Miocene climatic change might play an essential role in triggering these events. The onset of the Asian monsoon around the Oligocene/Miocene transition created a connection between forests from low to high latitudes of East Asia , which would provide a corridor for the MRCA of sections Thibeticae, Axiparis, and Euthyra to colonize subtropical East Asia. The enhancement of Asian summer monsoon since the middle Miocene established a humid climate in subtropical East Asia [84-86], and caused a significant expansion of forests in East Asia [83, 87]. These climatic and environmental shifts would create favorable habitats that facilitated the dispersal and divergence of sections Thibeticae, Axiparis and Euthyra in subtropical East Asia. Also, along with the expansion of forests in high latitudes of East Asia, the MRCA of P. quadrifolia and P. incompleta may have migrated into Europe.
The LTT analysis revealed that clade diversification within Paris abruptly accelerated since the late Miocene, which could be driven by the further strengthening of monsoonal climate in the summer and the initiation of the two distinct monsoon regimes that have gradually become established in subtropical East Asia since the late Miocene [83, 87, 88]. From then on, eastern, central and southern China and northern Indochina have been primarily governed by Pacific monsoon, whereas southwestern China and the Himalayas have been mainly affected by Indian monsoon [87–90]. The S-DIVA analysis (Fig.) and molecular dating (Fig. 3) showed that vicariance events occurred independently in sections Axiparis (9.36 Mya) and Euthyra (7.58 Mya) in the two regions mentioned above. This implies profound ecological heterogeneity resulting from climate differentiation may have driven significant allopatric speciation in the two regions [91–93]. In addition, the LTT analysis (Fig. 4) revealed that most extensive divergence in Paris, which was responsible for appearance of more than half of the extant taxa, took place in the Pliocene and the Pleistocene. It is believed that the Qinghai-Tibet Plateau (QTP) rose dramatically from the Late Miocene (ca. 10~8 Ma) to the early Pliocene (ca. 3.6 Ma) [94, 95], which dramatically modified global climate [89, 96] and thereby profoundly influenced biological processes in East Asia . During the Pleistocene, there were at least four major glaciations in East Asia , and these probably created significant isolation and diverse habitats in East Asia [99, 100]. Such complex geological, ecological, and environmental heterogeneity is expected to have driven diversification of a wide spectrum of plant clade [99, 101–104] and would also have facilitated a species radiative in Paris.
A negative correlation between genus-level diversity and the genus-average genome size was observed in plants [105, 106]. Knight et al.  proposed the large genome constraint hypothesis, which states that plant taxa with large size genomes diversify more slowly. Subsequently, Suda et al.  found that many island clades of Macaronesian angiosperms that underwent adaptive radiations have small genome sizes, and assumed that rapid diversification is more likely to happen in angiosperms with small genomes size. It is noteworthy that Paris is fairly distinctive in angiosperms in possessing large genomes. The minimum documented genome size in the genus (P. verticillata, 1C = 30.52 Gb) is much larger than the mean genome size (1C = 5.7 Gb) of angiosperms [108, 109]. Moreover, the known largest eukaryotic genome, that of P. japonica, 1 C = 148.88 Gb, belongs to Paris [76, 77]. In this study, we found that Paris may have undergone a species radiation since the late Miocene, which is not consistent with prediction that large genome size could limit speciation [106, 107]. It also has to be admitted that although there is sharp rise in the lineage diversification, the total number of species involved is not large in comparison to other radiations, for instance, Dianthus in the Mediterannean , and Aizoaceae in South Africa . The generality of the large genome constraint hypothesis needs to be further evaluated, although the increased lineage diversification detected here in Paris does not pose a major contradiction to it.