DOI: 10.5061/dryad.2jm63xsq5
title: Phylogenomic conflict analyses in the apple genus Malus s.l. reveal widespread hybridization and allopolyploidy driving diversification, with insights into the complex biogeographic history in the Northern Hemisphere

description: Phylogenomic evidence from an increasing number of studies has
 demonstrated that different data sets and analytical approaches often
 reconstruct strongly supported but conflicting relationships. In this
 study, 785 single-copy nuclear genes and 75 complete plastomes were used
 to infer the phylogenetic relationships and estimate the historical
 biogeography of the apple genus Malus sensu lato, an
 economically important lineage disjunctly distributed in the Northern
 Hemisphere and involved in known and suspected hybridization and
 allopolyploidy events. The nuclear phylogeny recovered the monophyly
 of Malus s.l. (including Docynia); however, the
 genus was supported to be biphyletic in the plastid phylogeny. An ancient
 chloroplast capture event in the Eocene in western North America best
 explains the cytonuclear discordance. Our conflict analysis demonstrated
 that ILS, hybridization, and allopolyploidy could explain the widespread
 nuclear gene tree discordance. One deep hybridization event (Malus
 doumeri) and one recent event (Malus coronaria) were detected
 in Malus s.l. Furthermore, our historical biogeographic
 analysis integrating living and fossil data supported a widespread East
 Asian-western North American origin of Malus s.l. in the
 Eocene, followed by several extinction and dispersal events in the
 Northern Hemisphere. We also propose a general workflow for assessing
 phylogenomic discordance and biogeographic analysis using deep genome
 skimming datasets.
description: Taxon
 sampling, library preparation, and deep genome skimming
 sequencing Taxon
 sampling is designed to resolve the phylogenetic placements and
 relationships of major clades within <em>Malus</em> s.l. For
 the convenience of discussion, we followed the widely accepted taxonomic
 system proposed by Phipps et al. (1990), in which they recognized five
 sections. Due to the great economic importance of
 <em>Malus</em>, numerous hybrid cultivars have been used as
 crops and ornamentals, which have hybridized either within one section or
 between sections, and this made it difficult to infer the phylogenetic
 relationships between cultivated and wild species. Hence, we excluded the
 widely recognized artificial hybrid species in <em>Malus</em>,
 e.g., <em>M. × asiatica</em> Nakai, <em>M. ×
 astracanica</em> hort. ex Dum.Cours., <em>M. ×
 cerasifera</em> Spach, <em>M. × dawsoniana</em> Rehder,
 <em>M.</em> <em>× domestica</em> (= <em>M.
 pumila</em> Mill.), <em>M.</em> <em>×
 floribunda</em> Siebold ex Van Houtte, <em>M.</em>
 <em>× halliana</em>, <em>M.</em> <em>×
 micromalus</em>, <em>M. × heterophylla</em> Spach,
 <em>M. × magdeburgensis</em> Schoch ex Rehder,
 <em>M.</em> <em>× prunifolia</em> (Willd.) Borkh.,
 <em>M. × sargentii</em> Rehder, <em>M. ×
 scheideckeri</em> Späth ex Zabel, <em>M. ×
 soulardii</em> (L.H.Bailey) Britton, <em>M.</em>
 <em>× spectabilis</em> (Aiton) Borkh., and <em>M. ×
 sublobata</em> Rehder. We hence sampled 39 ingroup individuals
 representing 18 wild species (out of ca. 24<sup>15</sup>),
 representing all five sections and <em>Docynia</em> within
 <em>Malus</em> s.l. In addition, due to the potential biphyly
 of <em>Malus</em> based on complete
 plastomes<sup>21,103</sup>, we sampled 38 outgroups across the
 tribes Maleae and Gillenieae to resolve the <em>Malus</em>
 phylogeny and identify possible events of cytonuclear discordance
 (Supplementary Table S1). We investigated 77 individuals in total, of
 which 27 species of deep genome skimming data were generated for this
 study (Supplementary Table S1). Total genomic DNAs were extracted
 from silica-gel dried leaves or herbarium specimens using a modified CTAB
 (mCTAB) method<sup>104</sup> in the lab of the Institute of
 Botany, Chinese Academy of Science (IBCAS) in China. The libraries were
 prepared in the lab of Novogene, Beijing, China using
 NEBNext<sup>®</sup> Ultra<sup>TM</sup> II DNA
 Library Prep Kit, and then paired-end reads of 2 × 150 bp were generated
 on the NovoSeq 6000 Sequencing System (Novogene, Beijing; 5-10 G data for
 each sample: Supplementary Table S1). Single-copy
 nuclear marker development The SCN marker development
 followed the pipeline in Liu et al. (2021). Briefly, the coding regions of
 <em>Malus domestica</em> (GenBank assembly accession:
 GCA_000148765.2) were first input into MarkerMiner
 v.1.0<sup>105</sup> to identify the putative single-copy
 genes. The resulting genes were then filtered by successively
 BLASTing<sup>106-108</sup> against six available genomes
 [<em>Malus baccata</em> (accession: GCA_006547085.1),
 <em>M. domestica</em>, <em>Pyrus betulifolia
 </em>Bunge (accession: GCA_007844245.1), <em>P. bretschneideri
 </em>Rehder (accession: GCA_000315295.1), <em>P.
 ussuriensis</em> Maxim. × <em>P. communis </em>L.
 (accession: GCA_008932095.1), and <em>P. pyrifolia
 </em>(Burm.f.) Nakai (accession: GCA_016587475.1)] in Geneious
 Prime<sup>109</sup>, with the parameters settings in the
 Megablast program<sup>110</sup> as a maximum of 60 hits, a
 maximum E-value of 1 × 10<sup>-10</sup>, a linear gap cost, a
 word size of 28, and scores of 1 for match and -2 for mismatch in
 alignments. We first excluded genes with mean coverage &gt; 1.1 for
 alignments, which generally indicate potential paralogy of the genes
 and/or the presence of highly repeated elements in the sequences. The
 remaining alignments were further visually examined to exclude those genes
 receiving multiple hits with long overlapping but different sequences
 during the BLAST. It should be noted that the alignments with mean
 coverage between 1.0 and 1.1 were generally caused by the presence of tiny
 pieces of flanking intron sequences in the alignments. These fragments
 were still accepted as a SCN gene here. After filtering, the remaining
 genes were used as references in the following gene assembly. The baits
 <em>in silico</em> could be available from the Dryad Digital
 Repository (Data 1): https://doi.org/10.5061/dryad.2jm63xsq5.
 Data
 processing and the assembly of single-copy nuclear genes Read processing and assembly
 followed the pipeline in Liu et al. (2021). Generally, we used Trimmomatic
 v. 0.39<sup>111</sup> for quality trimming and adapter
 clipping with the parameters (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:1:true
 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36). Then, the results were
 quality-checked with FastQC v. 0.11.9<sup>112</sup>. The
 number of clean reads after trimming was also calculated here for
 comparison (Supplementary Table S1). We used the HybPiper pipeline v.
 1.3.1<sup>113</sup> for targeting SCN genes with default
 settings; BWA v. 0.7.1<sup>114</sup> to align and distribute
 reads to target genes; SPAdes v. 3.15.0<sup>115</sup> with a
 coverage cutoff value of 8 to assemble reads to contigs; and Exonerate v.
 2.2.0<sup>116</sup> to align assembled contigs to target
 sequences and determine exon-intron boundaries. Python and R scripts
 included in the HybPiper pipeline<sup>113</sup> were used to
 retrieve the recovered gene sequences, summarize and visualize the
 recovery efficiency. Nuclear datasets
 construction and phylogenetic analysis Sequences for each SCN were
 aligned in MAFFT v. 7.480<sup>117</sup> with options
 “--localpair --maxiterate 1000”. Due to the variable sequencing coverage
 of each sample in this study, we employed three steps to remove the poorly
 aligned regions. We used trimAL v. 1.2<sup>118</sup> to trim
 the alignment of each SCN, in which all columns with gaps in more than 20%
 of the sequences or with a similarity score lower than 0.001 were removed.
 Given the low-quality assembly in some sequences,
 Spruceup<sup>119</sup> was used to discover, visualize, and
 remove outlier sequences in the concatenated multiple sequence alignments
 with the window size 50 and overlap 25. Because the Spruceup algorithm
 works better, the more data it has, we concatenated all the SCN gene
 alignments with AMAS v. 1.0<sup>120</sup> before running
 Spruceup. We also used AMAS v. 1.0<sup>120</sup> to split the
 processed/trimmed alignment back into single-locus alignments. The
 resulting alignments for each SCN were trimmed again using trimAL v.
 1.2<sup>118</sup> with the same parameters described above.
 Thirdly, we excluded the sequences less than 250 bp in each alignment with
 our customized python script (exclude_short_sequences.py, which can be
 available from Dryad Digital Repository
 https://doi.org/10.5061/dryad.2jm63xsq5) for decreasing the effect of
 missing data, because the short sequences in each alignment have few
 informative sites for the following coalescent-based species tree
 inference. The resulting SCN genes were used to infer individual ML gene
 trees using RAxML 8.2.12<sup>121</sup> with a GTRGAMMA model
 and the option “-f a” and 200 BS replicates to assess clade support for
 each SCN. TreeShrink v. 1.3.9<sup>122</sup> was used for
 detecting abnormally long branches in each tree with the default false
 positive error rate α = 0.05 and per-species mode. The shrunk trees and
 sequences have been used for the following phylogenetic inference, and
 hereafter these resulted sequences were referred to as “clean
 genes”. We
 generated three different datasets to reconstruct the phylogeny to account
 for the effect of missing data in each SCN gene: (1) 50%-sample dataset:
 each SCN gene with at least 900 bp and more than 50% samples (≥ 39
 individuals); (2) 80%-sample dataset: each SCN gene with at least 900 bp
 and more than 80% samples (≥ 62 individuals); (3) all-sample dataset: each
 SCN with at least 900 bp and more than 100% samples (77 individuals).
 These three datasets can be available from the Dryad Digital Repository
 (Data 2, 3, 4): https://doi.org/10.5061/dryad.2jm63xsq5. We used both
 concatenated and
 coalescent-based methods for phylogenetic inference of each dataset. We
 used PartitionFinder2<sup>123,124</sup> to estimate the
 best-fit partitioning schemes and/or nucleotide substitution models under
 the corrected Akaike information criterion (AICc) and linked branch
 lengths, as well as with rcluster<sup>125</sup> algorithm
 options for the nuclear dataset. The resulting partitioning schemes and
 evolutionary models were used for the following Maximum Likelihood (ML)
 tree using IQ-TREE2 v. 2.1.3<sup>126</sup> with 1000 SH-aLRT
 and the ultrafast bootstrap replicates and RAxML
 8.2.12<sup>121</sup> with GTRGAMMA model for each partition
 and clade support assessed with 200 rapid bootstrap (BS) replicates. The
 shrunk trees from TreeShrink<sup>122</sup> were used to
 estimate a coalescent-based species tree with ASTRAL-III (Zhang et al.,
 2018) using local posterior probabilities (LPP)<sup>127</sup>
 to assess clade support. Each of the gene trees was rooted, and low
 support branches (≤ 10) were collapsed using Newick
 Utilities<sup>128</sup> or
 <em>phyx</em><sup>129</sup> since collapsing gene
 tree nodes with BS support below a threshold value will help to improve
 accuracy<sup>130</sup>. In total, nine phylogenies were
 generated for topological comparison, and these nine trees are available
 from the Dryad Digital Repository:
 https://doi.org/10.5061/dryad.2jm63xsq5. Detecting and
 visualizing nuclear gene tree discordance To explore the discordance among
 gene trees, we employed <em>phyparts</em> v.
 0.0.1<sup>29</sup> to calculate the conflicting/concordant
 bipartitions by comparing the nuclear gene trees against the ML tree
 inferred from RAxML with a BS threshold of 50 (i.e., gene-tree
 branches/nodes with less than 50% BS were considered uninformative) for
 filtering out poorly supported branches, thus alleviating noise in the
 results of the conflict analysis<sup>29</sup>. We also used
 the internode certainty all (ICA) value that resulted from
 <em>phyparts</em> to quantify the degree of conflict on each
 node of a species tree given individual gene
 trees<sup>131</sup>. <em>Phyparts</em> results
 were visualized with <em>phypartspiecharts.py</em> (by Matt
 Johnson, available from https://github.com/mossmatters/MJPythonNotebooks/blob/master/phypartspiecharts.py). Furthermore, in order to distinguish lack of support from conflicting support in the species tree, we conducted Quartet Sampling (QS)<sup>41</sup> analysis with 100 replicates and the log‐likelihood cutoff 2. The QS method subsamples quartets from the input tree and alignment to assess the confidence, consistency, and informativeness of internal tree relationships, and the reliability of each terminal branch, and then four values are given in this analysis: QC = Quartet Concordance, QD = Quartet Differential, QI = Quartet Informativeness, and QF = Quartet Fidelity. The QS result was visualized with plot_QC_ggtree.R (by ShuiyinLIU, available from https://github.com/ShuiyinLIU/QS_visualization). Both the <em>phyparts</em> and QS results can provide alternative evidence for evaluating the discordance between gene trees. Comparing the nine topologies inferred above, we found two species (<em>Malus coronaria</em> and <em>M. kansuensis</em>) with conflicting phylogenetic positions among trees (referring to the result below). We used the “alternative relationship test” in<em> phyckle</em><sup>28</sup> to investigate conflicting bipartitions and discover the gene trees supporting each bipartition. Considering the complex origin of <em>Malus coronaria</em>, we employed the two bipartitions supported by Phylonetworks analysis (see Fig. 4 & Table 1 in Results). Two or more user-specified alternative bipartitions could be used as a constraint to infer gene trees. Arbitrarily, we set the cutoff of ∆lnL &gt; 100 as the outlier genes<sup>132</sup>. The resulting gene dataset supporting each bipartition was then used to estimate phylogenetic inference based on the concatenated (IQ-TREE2 and RAxML) and coalescent-based methods (ASTRAL-III) mentioned above. The resulting 12 trees could be available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.2jm63xsq5. Coalescence simulation and phylogenetic network estimation To measure the goodness-of-fit of the coalescent model with ILS explaining the gene tree discordance adequately, we performed a coalescent simulation analysis following the methods in previous studies<sup>34,133-135</sup>. Briefly, if the simulated gene trees based on the coalescent model correspond well to the empirical gene trees, the gene tree discordance may be explained by ILS. Given the calculation of gene tree distances, we subsampled 27 species representing each major clade in the phylogenetic framework of the tribe Maleae. The function “sim.coaltree.sp” implemented in the R package Phybase v. 1.5<sup>136</sup> has been used to simulate 10,000 gene trees with the multispecies coalescent (MSC) model. Hybrid Detector (HyDe) can be used to detect hybridization using phylogenetic invariants arising under the coalescent model with hybridization<sup>137</sup>. We sampled 40 taxa, including all 39 <em>Malus</em> s.l. individuals and one outgroup (<em>Pyrus communis</em>) to detect the possible hybridization events within <em>Malus</em> s.l. This dataset can be available from the Dryad Digital Repository (Data 6): https://doi.org/10.5061/dryad.2jm63xsq5. The γ denotes the inheritance probability of parent 1 (P1), while 1 - γ would be the probability of the hybrid population being sister to parent 2 (P2). Generally, significant γ values close to 0.5 indicate a recent hybridization event; significant γ values closer to 0 or 1 indicate an ancient hybridization event remained in the extant species. We herein set the γ threshold at 0.3 and 0.7, which followed convention<sup>132</sup>. To explore the possibility of reticulation as a cause of discordance in the apples and their allies, we employed the Species Networks applying Quartets (SNaQ) method<sup>43</sup> implemented in the software PhyloNetworks 0.14.0<sup>44</sup>, which explicitly accommodates introgression/gene flow and ILS. Given the computational limitation of PhyloNetworks, we used two datasets to test for hybridization events, i.e., 27-taxa sampling at the tribe level of Maleae and 14-taxa sampling at the genus level of <em>Malus</em>. These two datasets are available from the Dryad Digital Repository (Data 7 & 8): https://doi.org/10.5061/dryad.2jm63xsq5. Considering that <em>Malus</em> member<em>s</em> have been reported to have hybridized with many genera in Maleae, e.g., <em>Aria </em>(Pers.) Host and<em> Torminaria</em> M.Roem.<sup>5</sup>, we sampled 27 individuals, including five species representing each major clade in <em>Malus</em> s.l. and 22 outgroup species in Maleae. This taxon sampling scheme represents a reasonable compromise between taxonomic coverage and computational cost. The 27-taxa dataset construction followed the method described above (i.e., Nuclear dataset construction and phylogenetic analysis). The best trees generated from RAxML were used to estimate the quartet concordance factors (CFs), representing the proportion of genes supporting each possible relationship between each set of four species. The resulting CFs and the ASTRAL species tree were used as initial input to run SNaQ analysis (<em>h </em>= 0), and the resulting best network was used as starting topology to run the next <em>h</em> value (<em>h </em>+ 1), and so on. We investigated <em>h</em> values ranging from 0 to 5 with 50 runs in each <em>h</em> for estimating the best phylogenetic network. Each run generated a pseudo-deviance score: a value for fitting the network to the data, and estimated the inheritance probabilities (i.e. the proportion of genes contributed by each parental population to a hybrid taxon) for each network. Similarly, we also sampled 14 species, including 13 <em>Malus</em> species and one outgroup (<em>Pyrus communis</em>), to test the hybridization events among <em>Malus</em> members. The method followed that of the 27-taxa sampling dataset mentioned above. The best network was visualized using Dendroscope v 3.7.4<sup>138</sup>. Plastome assembly, annotation, phylogenetic analysis, and cytonuclear discordance A two-step strategy was used for obtaining high-quality chloroplast genomes. NOVOPlasty v. 4.3.1<sup>139</sup> was applied first to assemble the plastomes with high-quality raw data. Then we used the successive assembly approach<sup>140</sup>, combining the reference-based and the <em>de novo</em> assembly methods to assemble the remaining low-quality samples. With the <em>de novo</em> assembly and a seed-and-extend algorithm, NOVOPlasty was the least laborious approach and resulted in accurate plastomes; however, this program needs sufficient high-quality raw reads without gaps to cover the whole plastome. The whole plastomes assembled from NOVOPlasty then could be used as references for assembling the remaining samples. The successive method provided an excellent approach to obtaining relatively accurate and nearly complete plastomes with or without gaps from lower-coverage raw data. Due to the sensitivity of Bowtie2 v. 2.4.2<sup>141</sup> to the reference, this successive method needs a closely related reference sequence with increased time and RAM requirements. Several recent studies have described the procedure in detail<sup>21,60,103,140,142,143</sup>. All assembled plastomes have been submitted to GenBank with the accession numbers listed in Supplementary Table S1. The assembled plastid genomes from the low-coverage and high-coverage datasets were annotated using PGA<sup>144</sup> with a closely related plastome (MN062004: <em>Malus ioensis</em>) downloaded from GenBank as the reference, and the results of automated annotation were checked manually. The coding sequences of plastomes were translated into proteins to manually check the start and stop codons in Geneious Prime<sup>109</sup>. The custom annotations in the GenBank format were converted into the FASTA format, and five-column feature tables file required by NCBI submission using GB2sequin<sup>145</sup>. Given the considerable variation among plastid introns at the level of Maleae, we extracted 80 coding genes (CDSs) using Geneious Prime<sup>109</sup>, and these CDSs were aligned by MAFFT v. 7.475<sup>117</sup> with default parameters, respectively. This dataset with 77 samples and 80 plastid coding genes is available from the Dryad Digital Repository (Data 5): https://doi.org/10.5061/dryad.2jm63xsq5. The best-fit partitioning schemes and/or nucleotide substitution models for each dataset were estimated using PartitionFinder2<sup>123,124</sup>, under the corrected Akaike information criterion (AICc) and linked branch lengths, as well as with rcluster<sup>125</sup> algorithm options. The partitioning schemes and evolutionary model for each subset were used for the downstream phylogenetic analysis. Like the nuclear analysis, we estimated the ML tree by IQ-TREE2 v. 2.1.3<sup>126</sup> with 1000 SH-aLRT and the ultrafast bootstrap replicates and RAxML 8.2.12<sup>121</sup> with GTRGAMMA model for each partition and clade support assessed with 200 rapid BS replicates. We also used ASTRAL-III<sup>130</sup> for estimating a coalescent-based species tree. These three trees are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.2jm63xsq5. Cytonuclear discordance at various levels was detected in <em>Malus</em> s.l. (Fig. 7); coalescent simulation was also employed for evaluating the importance of ILS in explaining cytonuclear discordance<sup>35,135</sup>. We used the R package Phybase v. 1.5 to simulate 10,000 gene trees based on the ASTRAL species tree from SCN genes under MSC model. We used <em>Phyparts</em> v. 0.01<sup>29</sup> to explore conflicts among the simulated gene trees and plastid tree, and the proportion of gene tree concordance compared to the plastid genome tree was used to qualify the discordance. Dating and ancestral area reconstruction We aim to estimate the age of divergence of the three major clades identified from hundreds of SCN genes and 80 plastid coding genes in <em>Malus</em> s.l. The fossils used in this study are listed in Table 3. <em>Malus obensis</em> and the clade <em>Malus doumeri</em> - <em>Docynia delavayi</em> were grouped together because of the fruit similarity. <em>Malus parahupehensis</em> was thought to be similar to the living species <em>M. hupehensis</em>; however, this fossil species showed more similarities to <em>M. doumeri</em> because of the dentate margin and parallel, craspedodromous venation<sup>15</sup>. Thus, we grouped <em>M. parahupehensis</em>, <em>M. doumeri,</em> <em>M. obensis</em>, and <em>Docynia delavayi</em> as monophyletic. The earliest fossil, <em>Malus kingiensis,</em> was found in the middle Eocene from the western Kamchatka peninsula and can not be assigned to any living lineages of <em>Malus</em> s.l. based on the morphology; thus, this fossil species likely represented the stem clade of <em>Malus doumeri</em> and <em>Docynia delavayi</em>. However, only one fossil, <em>Malus antiqua</em>, was found in Europe, and this fossil with deeply lobed leaves was similar to the Mediterranean species. Therefore, we grouped this fossil with the two Mediterranean living species (<em>M. florentina</em> and <em>Eriolobus trilobatus</em>). Five fossil species were described from North America, especially the abundant fossil records in western North America. More or less deeply lobed leaves characterized all these leaf fossils, significantly distinct from the living western North American species (<em>Malus fusca</em>). They showed more similarities to the eastern North American and Mediterranean species, and these species were thus grouped together. Additionally, leaf fossil of <em>Malus</em> or <em>Pyrus</em> from the Republic site, Washington<sup>146</sup> was used to constrain the divergence between <em>Malus</em> and <em>Pyrus</em> at 46-44 Mya. We ran the dating analyses based on the SCN and plastid datasets to test the divergence time differences with significant cytonuclear discordance. Given the intensive computational burden of dating analysis using BEAST2, we employed a 19-taxa dataset with only one individual for each species in <em>Malus</em> s.l. and <em>Pyrus communis</em> as the outgroup, and this dataset is available from the Dryad Digital Repository (Data 13): https://doi.org/10.5061/dryad.2jm63xsq5. The divergence time estimation was run under a GTR model with a gamma rate inferred from PartitionFinder2<sup>123,124</sup>, an uncorrelated lognormal relaxed clock<sup>147</sup>, and the fossilized birth-death model<sup>148,149</sup>. Markow Chain Monte Carlo (MCMC) chains were run for 100,000,000, sampling every 20,000 generations in five parallel jobs. We used the LogCombiner v1.10 to combine log and tree files from the five independent runs of BEAST. The MCMC trace file was analyzed in Tracer v1.7.1<sup>150</sup>. Maximum credibility trees were generated in TreeAnnotator v1.10, and FigTree v1.4.4 visualized the MCC tree. To test the ancestral areas of three major clades of <em>Malus</em>, we conducted the ancestral area construction using BioGeoBEARS v. 1.1.1<sup>151</sup> implemented in RASP v. 4.2<sup>152</sup>. Geological evidence suggests that an aridity barrier existed from the western-most part of China to the eastern Asian coast from the Paleogene to the Miocene. It has been hypothesized to have acted as a climate barrier between these two regions<sup>153</sup>; thus, we subdivided East Asia into the northern and the southern areas<sup>154-159</sup>. Six biogeographic areas were defined across the distribution of <em>Malus</em> s.l.: (A), Southern East Asia; (B), Northern East Asia; (C), Europe and Central Asia; (D), Mediterranean; (E), eastern North America; (F), western North America. The MCC tree summarized by TreeAnnotator was used as input of RASP. The <em>maxarea</em> was set to six, i.e., the number of potential areas of a hypothetical ancestor was restricted to a maximum of six regions. The model with the highest AICc_wt value has been chosen as the best model.