Chloroplast phylogenomic analysis provides insights into the evolution of the largest eukaryotic genome holder, Paris japonica (Melanthiaceae)

Background Robust phylogenies for species with giant genomes and closely related taxa can build evolutionary frameworks for investigating the origin and evolution of these genomic gigantisms. Paris japonica (Melanthiaceae) has the largest genome that has been confirmed in eukaryotes to date; however, its phylogenetic position remains unresolved. As a result, the evolutionary history of the genomic gigantisms in P. japonica remains poorly understood. Results We used next-generation sequencing to generate complete plastomes of P. japonica, P. verticillata, Trillium govanianum, Ypsilandra thibetica and Y. yunnanensis. Together with published plastomes, the infra-familial relationships in Melanthiaceae and infra-generic phylogeny in Paris were investigated, and their divergence times were calculated. The results indicated that the expansion of the ancestral genome of extant Paris and Trillium occurred approximately from 59.16 Mya to 38.21 Mya. The sister relationship between P. japonica and the section Euthyra was recovered, and they diverged around the transition of the Oligocene/Miocene (20 Mya), when the Japan Islands were separated from the continent of Asia. Conclusions The genome size expansion in the most recent common ancestor for Paris and Trillium was most possibly a gradual process that lasted for approximately 20 million years. The divergence of P. japonica (section Kinugasa) and other taxa with thick rhizome may have been triggered by the isolation of the Japan Islands from the continent of Asia. This long-term separation, since the Oligocene/Miocene boundary, would have played an important role in the formation and evolution of the genomic gigantism in P. japonica. Moreover, our results support the taxonomic treatment of Paris as a genus rather than dividing it into three genera, but do not support the recognition of T. govanianum as the separate genus Trillidium. Electronic supplementary material The online version of this article (10.1186/s12870-019-1879-7) contains supplementary material, which is available to authorized users.


Background
Angiosperms exhibit extreme diversity in genome size that is defined as the haploid nuclear DNA amount, varying by approximately 2400-fold between the smallest and largest genomes [1][2][3][4]. Although the distribution of genome size in angiosperms is strongly skewed towards small genomes (with a mean value of 1C = 5.7 Gb and a modal value of 1C = 0.6 Gb) [4], to date, five species with the genome size 1C > 100 Gb have been documented. These plant species belong to the monocotyledonous family Melanthiaceae (one species in Paris, two species in Trillium), Liliaceae (one species in Fritillaria), and eudicot family Viscaceae (one species in Viscum) [5][6][7][8][9], suggesting that genomic gigantism may have originated and evolved independently in only a few lineages [1,10].
Because of the technical challenges in sequencing very large or very small genomes, insights into the mechanisms that drive the formation of genomic gigantism remain limited [9]. High-resolution and well-supported phylogenetic relationships between species with giant genomes and their closely related taxa can build evolutionary frameworks to elucidate the evolutionary history of these genomic gigantisms [9][10][11][12]. Unfortunately, a robust phylogeny for the genera Paris, Trillium and Viscum, which include genomic gigantisms, remains elusive [13][14][15], which impedes our understanding of the mechanisms underlying the formation and evolution of giant genomes.
Although the genome size of Polychaos dubia, a unicellular eukaryote, has been estimated to be over 670 Gb [16], this measurement is considered unreliable and inaccurate [4]. To date, the confirmed largest genome in eukaryotes has been observed in Paris japonica (Franch. et Sav.) Franch. (also known as Kinugasa japonica (Franch. et Sav.) Tatew.et Sutô.), with the 1C value of 148.88 Gb [1,17]. This plant is a perennial herb belonging to the monocotyledonous family Melanthiaceae tribe Parideae [18,19], and occurs natively in central and northern Honshu, Japan [20,21]. Cytological studies revealed that P. japonica is an octoploid with a chromosome number of 2n = 8x = 40 [1,22,23]. Because of its distinctive characters, such as showy and white sepals, and octoploid chromosome count, P. japonica has been historically placed either in the genus Paris (section Kinugasa) [21,23] or treated as a monotypic genus Kinugasa [20,24]. Moreover, the evolutionary relationships of P. japonica with related taxa have remained controversial in recent analyses based on single or multiple-locus DNA sequences. An analysis using the plastid rbcL region indicated that P. japonica is a sister to the genus Trillium [25]. A combination analyses of the plastid rbcL and matK and nuclear ITS DNA regions revealed that P. japonica is closely related to the genus Daiswa (=Paris section Euthyra) [13,26]. By contrast, two independent studies that based on the plastid psbA-trnH, trnL-F and nuclear ITS sequence data [27], and the combination of five plastid regions (atpB, rbcL, matK, ndhF and trnL-F) [28], resolved P. japonica as the sister group of the section Paris. These conflicts suggest that the relationships between P. japonica and allied taxa require further investigation.
Phylogenetic analysis using too few DNA sequences may result in a conflict between different sequence regions [29,30]; in such a case, it is not possible to reconstruct a robust and reliable phylogeny, in particular, at low taxonomic levels [31]. Because of its high level of intra-and infra-specific sequence variation, complete plastome DNA sequeces can offer valuable information for the analysis of complex evolutionary relationships in plants [32][33][34]. With the advent of next-generation DNA sequencing technologies, plasotmes have been widely used in recent years to reconstruct robust phylogenies for several phylogenetically difficult plant taxa [31,[35][36][37]; these cases suggest that whole plastome sequencing may provide novel evidence to elucidate the relationships between P. japonica and allied taxa. Despite the fact that the analysis of maternally inherited DNA loci may not demonstrate the complete history of the species, the complete plastome-based phylogeny can give us some valuable information to elucidate the maternal origin and evolution of the genomic gigantisms in P. japonica.
In the current study, we used low-coverage genome shotgun sequencing [38] to generate plastomes of P. japonica, P. verticillata, Trillium govanianum, Ypsilandra thibetica and Y. yunnanensis and then inferred the molecular evolution by comparing the structure and gene content to those of other published plastomes in Melanthiaceae. Then, we reconstructed the evolutionary relationships within the family to investigate the phylogenetic position of P. japonica. Finally, we dated the divergence of P. japonica to provide insights into the evolutionary history of the largest eukaryotic genome holder.
Although the gene content and arrangement were almost identical, pseudogenization and gene loss were found to have occasionally occurred within the family Melanthiaceae. Because of the presence of several internal stop codons in coding regions, cemA was identified as a pseudogene in all Paris and Trillium plastomes (Fig. 1a). In addition, the loss of the first exon of rps16 gene was found in the plastomes of Veratrum patulum and Chionographis japonica (Fig. 1a). Expansion of the IR regions into the ycf1 gene at the IR/SSC boundary occurred identically in all plastomes in Melanthiaceae, whereas their IR/LSC junctions were significantly variable. Three types of IR/LSC boundaries were observed in Melanthiaceae and outgroup taxa (Fig. 1b). The expansion of IR into the trnH-rps19 intergenic spacer (type III) was only found in V. patulum, whereas the expansion of IR into rps19 (type II) occurred in Trillium cuneatum, T. maculatum, and Paris polyphylla var. chinensis, as well as in outgroup taxa. Comparatively, characterized by the IR/LSC boundary falling into rps3, type I was observed in the remaining taxa (Fig. 1a).
The length of the intergenic region between rpl23 and ycf2 exhibited substantial variation among plastomes in the family Melanthiaceae, within which single-copy, duplicates and triplicates of trnI-CAU were observed (Fig.  1c). Triplication of trnI-CAU (type C) was observed in P. quadrifolia and P. verticillata (section Paris), whereas duplication of trnI-CAU (type B) was found in T. maculatum. A single-copy of trnI-CAU (type A) was identified in the other plastomes (Fig. 1a).

Phylogenomic analysis and divergence estimation
The tree topologies from both ML and BI analyses were identical. The phylogenetic relationships among the plastomes are presented in Fig. 1a. Five well-supported clades (BS = 100%, PP =1), corresponding to the five tribes (Melanthieae, Chionographideae, Heloniadeae, Xerophylleae, and Parideae) recognized by Zomlefer [18], were recovered. The tribe Melanthieae was sister to the rest of Melanthiaceae (BS = 100%, PP =1). The sister relationships between Chionographideae and Heloniadeae, as well as between Xerophylleae and Parideae, were fully supported (BS = 100%, PP =1). The intra-tribe relationships from our phylogenomic analysis are congruent with those of previous studies based on the nuclear ribosomal ITS and plastid trnL-trnF regions [18]; the combination of plastid DNA sequences [28,39]; and the plastid genome sequencing [15].
Within the tribe Parideae, the sister relationship between Trillium and Paris was recovered (BS = 100%, PP =1). The Paris species were further grouped into three fully supported lineages (BS = 100%, PP =1) that correspond to either the three narrowly-defined genera (Paris s.s., Kinugas and Daiswa, respectively) by Takhtajan [24] or the three sections (section Paris, section Kinugasa and section Euthyra, respectively) circumscribed by Hara [23]. Among them, P. japonica (section Kinugasa) was sister to the section Euthyra (BS = 100%, PP =1), and the section Paris was sister to the clade consisting of section Kinugasa and section Euthyra. The intersectional relationships obtained here are consistent with those of a previous study [40].
Three calibration points in Melanthiaceae (Fig. 1a) suggested by previous study [41] were used to constrain the plastome-based phylogenetic tree. The results suggested that the most recent common ancestor (MRCA) for the tribe Parideae dated at approximately 59.

Discussion
Robust phylogenies for species with giant genomes and allied taxa can build evolutionary frameworks to elucidate the origin and evolution of these genomic gigantisms [9][10][11][12]. Previous studies [13,[25][26][27] revealed that it is difficult to reconstruct high-resolution and well-  supported phylogenetic relationships between P. japonica, the largest eukaryotic genome holder, and its allied taxa based on too few DNA sequence regions. In this study, we sequenced the whole plastomes of P. japonica, as well as P. verticillata, Trillium govanianum, Ypsilandra thibetica and Y. yunnanensis. Coupled with publicly available plastomes in Melanthiaceae, we performed comparative and phylogenetic analyses of whole plastomes to clarify the evolutionary relationships of P. japonica with its closely related taxa. This study gives us some new information about the origin and evolution of the genomic gigantisms in P. japonica.

Plastome comparison
The loss of the first exon of rps16 was observed in the phylogenetically distinctive tribes Melanthieae and Chionographideae. Furthermore, the loss of trnD-GUC was only found in Y. thibitica. These results support the deduction that the loss of certain plastid genes may have independently occurred over the evolutionary history of angiosperms [32,42]. Therefore, the loss of certain plastid genes may not provide relevant evolutionary information. However, neither gene loss nor gene relocation were observed in any of the Melanthiaceae plastomes, implying the gene content and plastome structure in the family are highly conserved.
Previous studies have revealed that the protein-coding gene cemA has been lost in several non-photosynthetic parasitic plants [43][44][45]. To our knowledge, pseudogenization of this gene in photosynthetic autotrophic angiosperms has been only detected in the closely related genera Paris and Trillium (Fig. 1a). Although its function remains unclear [46], this mutation may provide a molecular synapomorphy to recognize the tribe Parideae [47]. In addition, as proposed in a previous study [15], the lineage-specific triplication of trnI-CAU in P. quadrifolia and P. verticillata could be used as a molecular synapomorphy to circumscribe the section Paris (Fig. 1a).
The IR/LSC boundaries of monocot plastoms generally expand into the trnH-rps19 gene cluster and the IR expansion duplicate trnH gene, which differs from those of non-monocot angiosperms [47]. In this study, we identified three types of IR/LSC expansions within Melanthiaceae; of those, type II and III exhibited the typical monocot IR/LSC junctions, whereas the IR/ LSC junctions of type I fell in rps3. Although IR/LSC expansions into the rps19-rpl22 intergenic spacer or rpl22 have been observed in some monocot orders, such as Asparagales, Commelinales, Zinbiberales and Poales [48][49][50], the more progressive expansion of IR/ LSC into rps3 has only been found in Melanthiaceae to date. The phylogenetic distribution of the three types of IR/LSC boundary in the tree topology suggests that the type III can be the ancestral state in Melanthiaceae, by compared with the expansion of IR regions into rps3 occurring in the derived tribes such as Chionographideae, Heloniadeae, Xerophylleae, and Parideae (Fig. 1a). Furthermore, the observation of type II of IR/LSC junction in T. cuneatum, T. maculatum, and Paris polyphylla var. chinensis may have been resulted from a secondary slippage of IR regions from rps3 to rps19.

Phylogeny inferences
Our phylogenomic analysis recovered five wellsupported lineages (BS = 100%, PP = 1) within Melanthiaceae, which correspond to the five tribes recognized by Zomlefer [18]. The evolutionary relationships recovered in this study are consistent with those of previous investigations [18,28,40,51] but with higher branch support (BS = 100%, PP = 1). The results further justify that whole plastid genome sequencing can improve the phylogenetic resolution in a certain lineage [33,34].
Our expanded sampling of the plastomes in Parideae provided an opportunity to reconstruct a robust intrageneric phylogeny in the tribe. The basal divergence in Parideae occurred approximately 38.21 Mya, forming two fully supported lineages (Paris and Trillium) in the tree topology (BS = 100%, PP = 1). The two genera share synapomorphies, including a single whorl of netveined leaves presenting at a stem apex, a stem apex bearing a solitary flower, and a chromosome base number n = 5 [16]. Within the clade Paris, the three sections (section Paris, section Kinugasa, and section Euthyra) outlined by Hara [23] as well as the three narrowly defined genera Paris s.s., Daiwa and Kinugasa by Takhtajan [24] were each recovered as monophyletic clades with strong support (BS = 100%, PP = 1) in both the ML and BI analyses. Given that species in the Paris clade share the morphological synapomorphies of flowers and leaves, 4-to 15-merous compared with the trimerous condition of Trillium [27], we correspondingly prefer to accept the taxonomic treatment of Paris as a single genus [21,23] rather than in three separated genera [24].
Since a previous study had not included the plastome of P. japonica in its phylogenetic analysis, its evolutionary relationships with other Paris species remained unresolved [15]. Both ML and BI analysis identically indicated that P. japonica (section Kinugasa) is a sister to the section Euthyra, which is congruent with the analyses of the plastid rbcL, matK and trnL-trnF regions [13,26,40]. However, the relationships recovered by our data largely differ from the results of combination analysis of plastid psbA-trnH and trnL-F and nuclear ITS sequences [27], and plastid atpB, rbcL, matK, ndhF and trnL-F regions [28]. It is noteworthy that, the well-supported sister relationship between P. japonica and the section Euthyra (BS = 100%, PP = 1) recovered in this study, can be also justified by the morphological synapomorphies that they share, such as a thick rhizome and angular ovary, in contrast to the long and slender rhizome and rounded ovary species of the section Paris (Fig. 2). In addition, the unusual morphological characteristics of the species (i.e., the showy, white sepals, and octoploid chromosome number) justify the taxonomic treatment of P. japonica as a distinctive section within the genus Paris by Hara [23].
Our data not only recovered the evolutionary backbone in Paris but also offered evidence to clarify disputes about the phylogenetic position of T. govanianum, which occurs natively in the Himalayan mountains. Although T. govanianum has a trimerous flower and leaves like those of Trillium species, it shares morphological features, such as narrow sepals and filiform petals, with Paris species (Fig. 3). Accordingly, T. govanianum was recognized as a separate genus Trillidium [13,52]. However, neither the ML nor BI tree topology separated T. govanianum from the Trillium species but grouped them into a wellsupported clade (BS = 100%, PP = 1). It is notable that similar finding has been shown in the phylogenetic analysis based on five plastid DNA regions that has a more extensive taxon sampling of Melanthiaceae [28]. Taken together, the results suggest that T. govanianum should remain in the genus Trillium and deny the recognition of the genus Trillidium.

Insights into the origin and evolution of the genomic gigantism in Paris japonica
The robust phylogeny reconstructed in the current study provided insights into the origin and evolution of the genomic gigantism in P. japonica. Most species in Melanthiaceae possess small or very small genomes, while large or giant genomes have been exclusively found in the two genera: Paris and Trillium [40]. Character reconstruction revealed that a genome size increase (more than four-fold) possibly occurred after the divergence of Xerophylleae and Parideae, but before the differentiation between Paris and Trillium [40]. Molecular dating indicated that the stem age and crown age of Parideae were approximately 59.16 Mya and 38.21 Mya, respectively, suggesting that the massive genome expansion would have lasted for a long period of approximately 20 million years. During this period, the ancestral genome of extant Paris and Trillium would have gradually expanded, implying that the genome size increase in Parideae could be the slow accumulation over tens of millions of years as a previous study proposed [40].
The phylogenomic analyses indicated that the section Paris is sister to the clade including P. japonica (section Kinugasa) and the section Euthyra. The relationships suggest that the formation of a giant genome in P. japonica most likely took place after the divergence of the sections Euthyra and P. japonica. Except for P. japonica, two species (T. × hagae and T. rhombilolium) with genome sizes 1C > 100 Gb have been found in the genus Trillium [5,6,9]. As we did not obtain samples of these two plants, their phylogenetic positions within Trillium remain unclear. Nevertheless, the evolutionary relationships of P. japonica with related taxa recovered in the study reveal that the formation of the giant genomes in P. japonica and Trillium species may have been independent events.
The coalescence of the plastomes of P. japonica and the section Euthyra occurred around the transition of the Oligocene/Miocene (20.30 Mya, 95% HPD: 34.64-9.96 Mya), when the opening of the Japan Sea separated the Japan Islands from the continent of Asia [53]. Although P. japonica and the section Euthyra are closely related, they Fig. 2 Comparison of morphological features among Paris japonica (section Kinugasa), section Paris, and section Euthyra occupy distinct distributions: P. japonica is endemic to Japan, whereas species from the section Euthyra are chiefly distributed in subtropical China and the Himalayas [23]. Hence, the divergence of P. japonica and the section Euthyra may have been triggered by the isolation of the Japan Islands from the continent of Asia.
Notably, the genome size of P. japonica is approximately 2-3 folds larger than that of those species belonging to the section Euthyra [40]. A line of evidence justifies that the genome size variation in plants is under selective constrains and has not evolved by a pure drift process [54][55][56]. As a result, genome size can be strictly related to the environment and ecology of a species [57]. In general, plants with larger genomes share some morphological traits, such as large body and stomata size [58]. Due to the drought susceptibility of the plants with large stomata, only those species occurring in humid habitats can sustain larger genomes [56,58]. Compared with the monsoonal climate which is characterized by obvious precipitation seasonality in subtropical China and the Himalayas [59,60], the maritime climate of the Japan islands [61] would create relatively more humid habitats that facilitate the evolution of P. japonica toward genomic gigantism.

Conclusions
The evolutionary relationships of the largest eukaryotic genome holder, P. japonica, with its closely related taxa were investigated by comparative and phylogenetic analyses of their complete plastome DNA sequences. Comparative analysis across plastomes in Melanthiaceae revealed that their structures and gene contents are highly conserved and provided molecular synapomorphies for some lineages of Parideae. Phylogenomic analysis and molecular dating recovered the evolutionary backbone of Paris and thus elucidated the phylogenetic position of P. japonica. The tree topologies and molecular dating indicated that the expansion of the ancestral genome of extant Paris and Trillium was probably a gradual process lasting for approximately 20 million years; the divergence of P. japonica and the section Euthyra may have been triggered by the opening of the Japan Sea, which separated the Japan Islands from the continent of Asia around the transition of the Oligocene/Miocene (20.30 Mya). This long-term separation would have played an important role in the formation and evolution of genomic gigantism in P. japonica. The phylogenetic position of P. japonica implies that the giant genomes of Paris and Trillium may have formed and evolved independently, even though the two genera are closely related. In addition, our phylogenomic analysis strongly supports the taxonomic treatment of Paris as a genus rather than dividing it into three genera, but did not support the recognition of T. govanianum as the separate genus Trillidium.

Plant material and shotgun sequencing
Leaf tissues of P. japonica, P. verticillata, T. govanianum, Y. thibetica and Y. yunnanensis were collected in the field and then dried with silica gel (one individual per species). The vouchers were identified by Dr. Yunheng Ji and deposited at the herbarium of Kunming Institute of Botany, Chinese Academy of Sciences (KUN); the voucher information is presented in Table 3. Genomic DNA was extracted from~20 mg of leaf tissue using a modified CTAB method [62]. Approximately 5 μg of purified genomic DNA was sheared by sonication. Paired-end libraries with an average insert size 350 bp were prepared using a TruSeq DNA Sample Prep Kit (Illumina, Inc., USA) according to the manufacturer's protocol. Shotgun sequencing was performed on the Illumina HiSeq 2000 platform.

Plastome assembly, annotation and comparison
Raw Illumina reads were filtered by NGS QC tool kit [63] to remove adaptors and low-quality reads. The pipeline developed by Jin et al. [64] was used for de novo plastome assembly. The clean reads of Paris species, T. govaniamum and Ypsilandra species were mapped onto the reference plastomes of P. quadrifolia (Genbank accession: KX784051), T. tschonoskii (Genbank accession: KR780076) and Heloniopsis tubiflora (Genbank accession: KM078036) using the Bowtie v2.2.6 software [65] with its default parameters and preset options. All of the plastid-like reads were assembled into contigs by SPAdes v3.10.1 [66] with the k-mer defined as 75, 85, 95 and 105. A customized python script [64], which can use BLAST and a built-in library to search the plastid-like contig, was employed to connect verified contigs into plastomes in SPAdes v3.10.1 [66], with its default parameters. The results of de novo assembly were visualized and edited with Bandage v.8.0 [67].
The resulting plastomes were annotated by Dual Organellar Genome Annotator database [68]. The annotations were manually proofed using Geneious v10.2.3 [69]. The start and stop codons of protein-coding genes were checked manually. All of the identified tRNA was verified by tRNAscan-SE v1.21 [70], with the preset parameters. Functional classification of the plastid genes was determined by referring to the online database CpBase (http://rocaplab.ocean.washington.edu/old_website/tools/ cpbase). The maps of plastomes were constructed with the Organellar Genome DRAW program [71].
The general features of plastome, such as structural rearrangements, gene loss/pseudogenization, gene duplication, and expansion/contraction of the IR regions, have provided evolutionary information in previous studies [15,32,72]. Therefore, we performed comparisons of these features among Melanthiaceae plastomes. The gene content and arrangement were visualized and compared with the MUMmer 3.0 program [73]. The boundaries of the LSC, IR, and SSC regions in each plastome were compared using Geneious v10.2.3 [69].

Phylogenomic analysis
To examine the phylogenetic position of P. japonica, 24 plastomes representing wide phylogenetic diversity in the family Melanthiaceae were included in the phylogenomic analysis (Additional file 1: Table S1). The plastomes of Campynema lineare, Fritillaria cirrhosa, Luzuriaga radicans and Smilax china were used to root the tree. Of those, five plastomes were newly generated in the current study (Table 3), and the rest plastomes were downloaded from the NCBI database (Additional file 1: Table S).
The complete plastome DNA sequences were aligned using MAFFT [74] integrated in Geneious v.10.2.3 [69], with manual adjustment if necessary. The phylogenomic analyses were carried out with the standard Maximum Likelihood (ML) and Bayesian Inference (BI) methods. ML analyses were performed using RAxML-HPC Black-Box v8.1.24 [75] with 1000 replicates of rapid bootstrapping (BS) under the GTR-GAMMA model. The search of the best-scoring ML tree and rapid bootstrapping were performed in a single run. The choice of the best nucleotide sequence substitution model for BI analysis was determined using Modeltest v3.7 [51] with the Akaike Information Criterion [76]. BI was performed with MRBAYES v.3.1.2 [77] using the model (TVM + I + G) selected. Two independent parallel Markov Chain Monte Carlo (MCMC) runs with tree sampling every 100 generations for one million generations, with the first 25% discarded as burn-in, were conducted. Stationarity was considered to be reached when the average standard deviation of the split frequencies was < 0.01. The posterior probability values (PP) were determined from the remaining 0.75 million trees.

Molecular dating
To date, no fossils have been identified for the family Melanthiaceae and its close relatives. Calibrated by 17 fossils across the monocots and major clades of angiosperms, a previous study [41] revealed that the crown age of family Melanthiaceae was approximately 84.8 Mya, while the clades Parideae-Xerophyllideae and Chionographideae-Heloniadeae diverged approximately 74 Mya, and the tribes Parideae and Xerophyllideae split approximately 52.3 Mya. We used these events to calibrate the phylogenetic tree (Fig. 1a). Molecular dating was performed using the MCMCTREE v4.9c program integrated in the PAML program package [78]. The ML tree topology was used to estimate the divergence times of nodes. The independent-rates molecular clock was chosen as the clock model, and HKY85 was selected as the substitution model. The root age was set as less than 100 Mya. The divergence of Melanthiaceae was calibrated with a minimum age of 84.8 Mya. The node uniting Parideae-Xerophyllideae and Chionographideae-Heloniadeae was set to a minimum age of 74 Mya, while the divergence of Parideae and Xerophyllideae was set to a minimum age of 52.3 Mya. Other parameters were defined as their defaults. MCMC chains were run for 10,100,000 iterations. The first 100,000 iterations were discarded as burn-in, and trees were sampled every 10 iterations until 1000,000 samples were gathered.