Phylogeographic inference of Sumatran ranids bearing gastromyzophorous tadpoles with regard to the Pleistocene drainage systems of Sundaland

Rivers are known to act as biogeographic barriers in several strictly terrestrial taxa, while possibly serving as conduits of dispersal for freshwater-tolerant or -dependent species. However, the influence of river systems on genetic diversity depends on taxa-specific life history traits as well as other geographic factors. In amphibians, several studies have demonstrated that river systems have only minor influence on their divergence. Here, we assess the role of the paleodrainage systems of the Sunda region (with a focus on the island of Sumatra) in shaping the evolutionary history of two genera of frogs (Sumaterana and Wijayarana) whose tadpoles are highly dependent on cascading stream habitats. Our phylogenetic results show no clear association between the genetic diversification patterns of both anurans genera and the existence of paleodrainage systems. Time-calibrated phylogenies and biogeographical models suggest that these frogs colonized Sumatra and diversified on the island before the occurrence of the Pleistocene drainage systems. Both genera demonstrate phylogenetic structuring along a north–south geographic axis, the temporal dynamics of which coincide with the geological chronology of proto Sumatran and -Javan volcanic islands. Our results also highlight the chronic underestimation of Sumatran biodiversity and call for more intense sampling efforts on the island.

www.nature.com/scientificreports/ time periodically established connections between the different landmasses of the Sunda Shelf 18 and formed four extensive paleoriver systems in the Sunda and Sahul shelves, i.e., the Malacca Strait River System, the Siam River System, the North Sunda River System, and the East Sunda River System 19,20 . These Pleistocene drainage systems purportedly had impacts on the biodiversity patterns in the region as they could have served as potential dispersal routes between the Greater Sunda Islands 19,20 . For example, the freshwater riverine faunas of rivers presently restricted to Indo-China, the Malay Peninsula, or one of the Greater Sunda Islands were probably connected during the Pleistocene 11,19,21,22 . Extant species can be mapped onto the river systems and analyzed with biogeographic and/or phylogeographic approaches. The genetic structure of populations may have been strongly influenced by dispersal barriers resulting from fragmentation of the riverscape (defined as a mosaic of freshwater river habitat that is spatially structured and hierarchically organized across multiple scales 23 ). Genetic diversity may be partitioned and regionally distributed according to river drainages if the paleoriver systems played a significant role in the speciation and dispersal of the groups of organisms under study. This appears to be the case, at least for some freshwater organisms and other groups with strong aquatic dependencies. For example, de Bruyn et al. 12 found a strong correlation between phylogenetic structuring in clades of freshwater fishes and their distributions across Sundaland paleodrainages. Similar phylogeographic patterns have been observed in other taxa such as southern Indochinese amphibians 24 of the Lower Mekong in southern Indochina, and the Mekong mud snake 25 . The evidence from these studies suggests that the Quaternary landscapes of Indochina and the Sunda Shelf shaped the genetic divergence patterns in populations of certain taxa. The influence of rivers (either as corridors or barriers) on the distribution and genetic structure of local fauna has also been tested in various other taxa and regions (e.g., in the Amazon: frogs 26 , frogs and small mammals 27 , mammals 28 , birds 29 ; on Madagascar: mammals 30 , frogs 31 ; in the southern USA: fish 32 ; in eastern Australia: water skinks 33 ). These studies showed that the vicariant influence of river systems depends on taxon-specific life-history traits as well as the geographic setting 31 . Moreover, the current hydrography of river systems may not always explain observed species distribution patterns.
Most amphibians undergo a complex biphasic life cycle with larval forms that are restricted to and strongly depend on aquatic habitats, whereas the terrestrial stages are potentially more vagile and prone to dispersal. Still, rivers may serve as barriers in amphibian populations and may confine their dispersal [34][35][36] . Here we selected two endemic genera of Sumatran ranid frogs (Sumaterana and Wijayarana) to investigate whether the paleoriver systems of Sumatra played a role in structuring the distribution of their genetic diversity. Given the trans-island distribution of Wijayarana, we also included samples representing populations from the western part of Java. Species of Sumaterana and Wijayarana, the only Sumatran ranids that inhabit torrential streams, possess gastromyzophorous tadpoles [37][38][39] that are characterized by the possession of a large adhesive sucker at the abdomen as an adaptation to torrential stream habitats 40 . Because of their strong association with cascading streams during their larval stage, Sumaterana and Wijayarana species appear as suitable focal taxa to assess the role of paleodrainage systems of the Sunda region, in shaping the phylogenetic structure of regional stream-dependent anurans on the Sunda region. We particularly focused on Sumatra given that four major paleodrainage systems of Sundaland concurrently existed on this island 19 . Our results show that the distribution patterns of Sumaterana and Wijayarana do not correlate with the Quaternary paleodrainage systems of Sumatra. Instead, the diversification processes of these taxa appear to be more complex, with the focal taxa having colonized Sumatra and diversified therein much earlier than the formation of the Pleistocene River systems. More specifically, the spatiotemporal diversifications of both genera are best explained by the orogenesis of the Bukit Barisan mountain range and initial episodes of geographic isolation of landmasses (between the Miocene to Pliocene) that form the contemporary volcanic islands in the region.

Result
Phylogenetic relationships of Sumatran ranids with gastromyzophorous tadpoles. We generated a new sequence dataset (N = 146; sampling localities shown in Fig. S1), comprising ten concatenated genetic loci (mtDNA and nucDNA) for Sumatran frogs of the genera Sumaterana and Wijayarana along with sequences of closely related taxa within Ranidae (Supplementary Table S1). We used the final concatenated alignment of a total of 7,582 bp (28.91% proportion of missing information; see Supplementary Table S1) to infer phylogenetic relationships of Sumatran ranids that possess larvae with an abdominal sucker. Both maximum likelihood (ML) and Bayesian inference (BI) trees revealed the similar topologies for the relationships of ranid frogs with gastromyzophorous tadpoles (blue color in Fig. 1 Wijayarana sumatrana was closely related to W. sp3 from Java (BS/PP = 100/0.99; BS = 100). Minor differences between ML and BI trees also appeared in the intra-specific relationship of Meristogenys amoropalamus from Borneo. In the BI tree ( Supplementary Fig. S2), M. amoropalamus was sister to a clade of M. orphocnemis + M. poecilus, whereas in the ML tree this species was sister taxon of an unidentified Meristogenys species (Supplementary Fig. S3). Although the mainland Asian genus Amolops also possesses gastromyzophorous tadpoles, both trees suggest Amolops to be closely related to ranid species that do not possess these specialized tadpole types (e.g., Odorrana, Chalcorana, Pulchrana. Hylarana), rather than grouping with Huia, Meristogenys, Sumaterana, and Wijayarana that have this larval form. In both analyses (Fig. S2), Clinotarsus, a ranid taxon with non-gastromyzophorous tadpoles was recovered as nested within the clade of Sundaland ranids with gastromyzophorous tadpoles (Huia, Meristogenys, Sumaterana, Wijayarana), however, with low support values.
Sumatran paleodrainage systems and the distribution patterns of Sumatran ranids with gastromyzophorous tadpoles. We used the same BI tree (Fig. 1) to map the distribution of genetic variation The phylogeographic structure of Sumaterana and Wijayarana did not show any clear correlation to the paleodrainage systems on Sumatra (Fig. 2a, b; see "Methods" for paleodrainage systems definition). Members of both genera from each paleodrainage system (sensu Voris 19 ) east of the Barisan Mountain range were more closely related to members from their respective western watershed than to each other. Although Clade D shows a monophyletic group of samples from the Malacca Strait paleodrainage system, all samples were collected from different tributaries but within the same watershed (Sungai Rokan).
Considering their wide distribution across the island, we used Sumaterana crassiovis ( Fig. 2a Clade A-D), Wijayarana sp1 ( Fig. 2b Clade E-F), and W. sumatrana ( Fig. 2b Clade G-I) to hypothesize divergence scenarios in Sumatra. The three taxa are genetically structured into northern (Sumatra-North in Fig. 3) and southern (Sumatra-Central and Sumatra-South in Fig. 3) lineages. Clade A, E, and G are distributed in the northern part, whereas clades B2, C, F, H, and I exhibit southern distributions. The Sumatera Barat province (located in the center of the island and comprising parts of Sumatra-North and Sumatra-Central in Fig. 3) appears to be the zone of range overlap between northern and southern species distributions. Samples from the Sumatera Barat (SB) province occur in both groups. This is supported by the fact that clade B1, which comprises samples from the SB province, also contains samples from Jambi (JB) province. Clade I encompasses samples from Aceh (AC), Sumatera Barat (SB), and Lampung (LP) provinces.
Divergence time estimation and ancestral area reconstruction. We estimated divergence times of Sumatran frogs with gastromyzophorous tadpoles and their close ranid relatives to understand how the diversification history of these frogs relates to past geographic or climatic events on Sumatra. The chronogram (Fig. 3) resulting from this analysis differed slightly from the ML and BI trees in the arrangement of the terminal nodes (particularly for Wijayarana and Sumaterana). Furthermore, in the chronogram, Clinotarsus was placed as sister to a clade comprising taxa with gastromyzophorous tadpoles from Sundaland (Huia + Meristogenys + Sumaterana + Wijayarana) as opposed to being sister to Sumaterana in our ML and BI phylogenies.
Our inference of biogeographic history of both genera suggest that the clade comprising Huia + Meristogenys + Sumaterana + Wijayarana began to diversify in the Early Eocene (approximately 47.   www.nature.com/scientificreports/  Wijayarana have continued to diverge into many intra-specific lineages (Fig. 3).

Discussion
The biodiversity and diversification patterns of amphibians on Sumatra and most of the Sundaland in general, remain poorly understood 5,39,41 , despite a variety of phylogenetic and taxonomic studies in the last decade 39,[42][43][44][45][46] .
Here, we add to this existing knowledge by analyzing the evolutionary relationships and biogeography of two stream-adapted genera of Sumatran frogs using multi-locus genetic data and comprehensive sampling. In doing so, we corroborate the recognition of Wijayarana as a distinct genus with at least five independent lineages 39,47 . This discovery is consistent with the evidence that Sumatran amphibian diversity is gravely underestimated (see Figs. 2 and 3) and that cataloguing of species diversity in Sumatra would benefit greatly if more thorough sampling efforts were performed 39,[47][48][49] . Furthermore, the establishment of a stable taxonomy for Sumatran ranids with gastromyzophorous tadpoles (Sumaterana and Wijayarana), combined with the meticulous documentation  www.nature.com/scientificreports/ of their geographic distribution in Sumatra represent important steps towards a better understanding of regional diversification patterns and processes 50,51 .
The presence of the Pleistocene paleodrainage system in Sundaland has been proposed as one of the main factors driving diversification patterns in some taxa (e.g., Sundaland freshwater fish 12 and Mekong mud snakes 24 ). However, our results failed to garner evidence for paleodrainages sensu Voris 19 having significantly influenced evolution of Sumatran frogs Sumaterana and Wijayarana (Fig. 2), despite their strong dependence on rivers as larvae. In fact, the major diversification events in these frog lineages pre-date the formation of the paleodrainages by several million years (Fig. 3). Unlike previous studies, we also considered watershed systems that currently exist on the western slopes of the Barisan range in our biogeographic model (labeled with number 5 in Fig. 2). In this context, frogs from the eastern slopes of the Barisan Mountain range appear more closely related to individuals from the corresponding western slopes, despite belonging to distinct contemporary basins, which drain in opposite direction (i.e., west vs. east) (Fig. 2, Supplementary Table S1).
Unlike the freshwater fishes and the Mekong mud snakes that are restricted to fresh water, frogs with gastromyzophorous tadpoles presumably have two potential dispersal pathways: (1) via rivers during their larval stage and (2) via land during their adult stage. This implies that the dynamics of their spatio-temporal evolution are more complex than that of an organism with more restrictive means of dispersal (e.g., fish). Davis et al. 23 proposed that both distance and potential barriers to gene flow may shape the genetic variation of taxa. Building on this premise and based on our results, we hypothesized that the formation and configuration of the Barisan Mountain range (e.g., contour and topology) probably have had the greater influence as drivers of diversification in this group of frogs.
Overall, the diversification of frogs examined herein seems largely structured geographically along a north-south axis (i.e., the orientation of the Bukit Barisan mountain range) within Sumatra (Fig. 2). Members of the genus Sumaterana, for example, display a spatio-temporal diversification pattern that appears to advance from the north of the island towards its southern end. A similar phylogeographic scheme also presents itself in the results of de Bruyn et al. 3 who interpret this scenario as evidence for their "radiation by paleodrainage" diversification hypothesis for two Southeast Asian fish genera. However, the occurrence of this pattern of radiation in two vastly different local organisms (i.e., fish and frogs) at two vastly disparate temporal junctions (i.e., early Miocene vs. early Eocene onwards) indicates that the drivers of this phenomenon may be more numerous and more complex than just the configuration of Pleistocene River drainage systems.
The distribution boundary between the northern and southern groups appears to be located in the Sumatera Barat province (Fig. 2). However, some of the samples from this hypothesized boundary grouped either in northern or in southern lineages. Similar north-south genetic partitioning has also been reported for three species of Sumatran tree frogs of the genus Rhacophorus (R. catamitus, R. modestus, R. poecilonotus 13   www.nature.com/scientificreports/ influenced the population structure of these monkey species 54 and facilitated re-colonization on Sumatra 57 . However, O'Connell et al. 52 argued that the partition between the northern and southern populations of the Rhacophorus species was driven by the paleoclimate during the Miocene and Pleistocene epoch. A similar diversification process, with a slight variation, is also presented by Shaney et al. 58 for Sumatran montane dragons (genera Dendragama, Lophocalotes, and Pseudocalotes) wherein the lizards are genetically segregated into four major biographic zones: North, North-Central, Central, and South montane zones. Elevational boundaries, the paleoclimate during the Pleistocene, and ecological factors (e.g., competitive niche exclusion) have all been suggested as likely drivers of phylogenetic structuring in these agamid genera 58 . Thus, it remains hard to state with confidence whether the population genetic structuring shown by Sumaterana and Wijayarana (north-south partitioning) was solely influenced by the paleoclimate and geology of Sumatra (as in Rhacophorus 52 ), or also influenced by other factors (e.g., elevation, ecology as proposed by Shaney et al. 58 ). Given existing data, it is also unclear if the Toba eruption played a decisive role in the aforementioned spatio-temporal evolutionary dynamics of our focal taxa. Sumatra has had an extremely active and complex tectonic history and several local volcanic eruptions may have played a role the genetic partition and recolonization of regional species.
In the hypothetical scenarios of in-situ diversification of the Sumatran Rhacophorus, O'Connell et al. 52 suggested that the northern populations of these tree frogs were considered to have occupied the northern region of Sumatra at the beginning in the Miocene 52 when suitable habitats were likely to be abundant in the region 7 . According to our dated phylogeny (Fig. 3) and ancestral area reconstruction (Fig. 4) 4). These time estimations marked the beginning of in-situ diversification of the two genera on the island of Sumatra, which presumably occurred synchronously.
Hall 14,16 suggested that by the Mid-Late Oligocene major parts of Sundaland were subaerial including some freshwater lakes. During this time, the northern part of Sumatra had a more seasonal climate 8 and lacked active volcanoes 14, 16 . In contrast, the southern part of the island had a wet climate and harbored several active volcanoes. Consequently, habitat types in the northern and southern parts of Sumatra were very different. Although volcanoes were not present at the time in the northern region, some highlands, which provide ample cascading habitats in this area, were known to be present 16 . Our dated phylogeny and ancestral area reconstruction analyses suggest that the population of the most recent common ancestor (MRCA) of Sumaterana first evolved in cascading habitats of northern Sumatra (light brown color in the map; Fig. 3) after initial colonization from mainland at approximately 25.11 Ma. Considering some overlapping distributions of the three extant Sumaterana species, the MRCA of S. dabulescens and the MRCA of S. montana and S. crassiovis adapted to distinct but overlapping elevational ranges, thus occurring in sympatry in some areas (see Arifin et al. 39 ). In contrast, at ~ 27.65 Ma the MRCA of Wijayarana were migrated to Java from mainland Asia via Sumatra, subsequently experiencing an extinction event in its north and central ranges of Sumatra (see Fig. 4).
Between 20-15 Ma, the chain of volcanoes in Sumatra shifted northwest up to the northern part of Sumatra due to strike-slip faulting, and in response to the region-wide Sundaland deformation initiated by Australian plate's collision with eastern Indonesia 14, 16 . The increase of regional marine transgression resulted in increased sea levels and therefore large parts of Sumatra were gradually submerged 14,59 , leaving only high elevation mountains above water. Our ancestral area reconstruction suggests that the MRCA of Sumaterana dabulescens and the MRCA of S. montana and S. crassiovis (Fig. 3) could have been locally concentrated on different volcanic peaks serving as refugia on northern Sumatra and thus becoming geographically isolated from each other. The ever-wet climate across the whole island during this time 8 may have supported populations in limited numbers on mountain tops, impeding gene flow between isolated refugial populations. The same conditions may have applied to the MRCA of Wijayarana, which presumably occupied available refugia on mountain tops across the Miocene highlands of south Sumatra. On the other hand, the tectonic activity that shifted the position of the western-most end of Java to connect with the southernmost tip of Sumatra 14,16 , putatively allowed the MRCA of the Javan Wijayarana to colonize Sumatra at approximately 16.42 Ma onwards. This lineage eventually gave rise to W. sumatrana via vicariance due to the isolation enforced by a marine barrier. Around the same time, ancestral populations of Sumatran Wijayarana (Fig. 2) from the mountain-top refugia of southern Sumatra dispersed into the more central regions of the island.
As sea levels gradually decreased 16 from 10 Ma onwards, gene flow between previously geographically isolated populations could potentially have resumed. The MRCA of Sumaterana crassiovis could have then migrated further south along with the MRCA of S. montana, whereas the MRCA of S. dabulescens remained in the north. In contrast, the MRCA of Wijayarana sumatrana and the MRCA of W. sp1 and W. sp2 moved northwards occupying all suitable niches on the island. From 5 Ma (the beginning of the Pliocene) onwards, Sumatra experienced further land accretion, which dramatically modified the topology of the island. The number of volcanoes increased, leading to a greater number of high-altitude habitats in southern Sumatra 14, 16 . This geological process became more frequent during the Pleistocene and continued into the Holocene in many areas 60 . Seasonal climate occurred during this period 8 . Consequently, suitable cascading habitats for both Sumaterana and Wijayarana became more abundant across the island. Our overall results support the premise that the gradient of divergent environments between the northern and southern parts of Sumatra during the Pliocene onwards initiated the segregation of Sumatran fauna (e.g., Sumaterana and Wijayarana) into the northern and southern groups seen today (see Figs. 2 and 3). These environmental conditions could have been effective drivers of synchronous (insitu) diversification within the two genera. The diversification became more pronounced from the beginning of the Pliocene onwards, compared to the previous epoch (see Fig. 3).
In conclusion, our work is the first to leverage contemporary biogeographic approaches to elucidate the evolution of these Sumatran ranids with gastromyzophorous tadpoles. In doing so, our study sheds light on the www.nature.com/scientificreports/ diversification pattern of Sumaterana and Wijayarana on the island of Sumatra through space and time. The results of our study show that both Sumaterana and Wijayarana probably experienced rapid divergence from the Late Miocene or the Early Pliocene onwards (Fig. 3) driven primarily by the paleo-climatic (cooling and warming) and changing of land and sea level during this period 8,14,16 . However, we would also like to highlight that there are several caveats concerning this interpretation given the complex climatic and geological history as well as the lack of an adequate number of comparative studies for the region. Considering current distribution patterns for both genera, elevational gradients might have also been influential 58; unpubl. . Furthermore, ecological factors (e.g., competitive exclusion) may have also been drivers for the observed genetic segregation. For example, during our fieldwork, we often found Sumaterana and Wijayarana living in sympatry near the same stream. Nevertheless, Wijayarana tend to avoid the water whenever Sumaterana species were found in the stream pers. obs . Observation such as these indicate that the current evolutionary information on these ranids, would benefit greatly by having more complementary ecological data. Evolutionary and ecological investigations will both be vital for a clearer understanding of the complex biogeography of not only these unique amphibians but also of other Sundaland fauna in general.

Methods
Taxon sampling and molecular data. We sampled a total of 146 ranid frogs belonging to the genera    Supplementary Table S1. DNA extraction, PCR, and sequencing followed published protocols 39 . Sequences of ten loci were generated for five mtDNA (12S, 16S + tRNA val , COI, cyt b, ND2) and five nucDNA (Brain-derived neurotrophic factor, BDNF; Neurotrophin 3, NTF3; Proopiomelanocortin, POMC; recombination-activating gene 1, RAG1; tyrosinase exon 1, TYR) markers; see Supplementary Table S1 for Genbank accession numbers. Primer information and PCR annealing temperatures applied for this study are provided in Supplementary Table S2.
Phylogenetic analyses. The dataset of ten concatenated genes consisted of 7,582 bp. In order to create a dataset with a minimum amount of missing data, we only used individuals for which at least three loci (or 1,942 bp total length) were successfully sequenced (see Supplementary Table S1 for details on marker coverage for each sample). Our final dataset comprised 146 specimens for Sumaterana (N = 85) and Wijayarana (N = 61), and an additional 21 sequences for closely related taxa. We tested a variety of models and partitioning strategies to find the best partitioning scheme and substitution models for the concatenated dataset using PartitionFinder v.1.1 62 66 . In RAxML, eleven distinct partitions and associated models (GTR + Γ + G) were defined, and we performed joint branch length optimization. Tree support was obtained running 1,000 bootstrap replicates.
For the Bayesian analysis, we performed two independent runs with one cold and three heated chains for 50 million Markov Chain Monte Carlo (MCMC) iterations, sampling every 1,000 generations. Convergence of runs was assessed using the trace plot generated from MrBayes, the average split frequencies being < 0.01, and by assessing ESS values (> 200) of the log files with Tracer v.1.6 67  www.nature.com/scientificreports/ Mapping of distribution patterns. The BI tree was transformed and edited in MESQUITE v.3.6 71 into a color-coded tree to map the distribution patterns of the Sumatran ranid frogs with gastromyzophorous tadpoles. We followed Voris' 19 definition of the four paleodrainage systems in Sundaland (color-coded and labeled with number 1-4 in Fig. 2). Voris 19 did not specify any watersheds that ran into the Indian Ocean (western slope of the Bukit Barisan mountain range) in his paleodrainage definition. Our samples were also collected from both the western and eastern slopes of the Bukit Barisan mountain range. Thus, in this study, we categorized all watersheds on the western slopes of the mountains as distinct drainage system (color-coded and labeled with number 5 in Fig. 2).

Divergence time estimation.
We estimated divergence times of Sumatran ranids with gastromyzophorous tadpoles using BEAST v.2.4.8 72 . We performed analyses with the complete dataset (N = 167) and the partitioning scheme as described in Supplementary Table S3. The molecular clock was calibrated using two evolutionary rates from previous studies: (1) 1.41% (0.705% per lineage)/MY, a rate estimated for the cytb region of Sylvirana latouchii 73 , another ranid species; and (2) 1.00% (0.50% per lineage)/MY, as suggested by Kakehashi et al. 74 for the 16S region of anurans more generally. Though we used the same values for the mean and standard deviation as in Tominaga et al. 75 the prior distribution was set to a uniform distribution rather than a normal one (compared to Tominaga et al. 75 ) because the analyses failed to reach convergence under a normal distribution. First, we performed a series of test runs to find the optimal settings for our data. We varied the clock model, tree prior, and gamma hyperparameter and performed runs while retaining the remaining parameters as default. Each preliminary run was performed for 100 million iterations, sampling every 10,000 generations. We then compared the log files and chose the parameter set with the highest likelihood for our final analysis (see Supplementary Table S4 for likelihood comparison). The best parameter set was as follows: uncorrelated log-normal relaxed clock model, Yule tree prior, and HKY for the substitution model with four gamma categories. We performed the final analysis in two independent runs of Markov chains for 500 million generations, sampling every 10,000 generations. We used Tracer v. Ancestral area reconstruction. We performed statistical ancestral area reconstruction for the taxa of interest using BioGeoBEARS 76 as implemented in RASP 4.2 77 . We pruned our calibrated trees from BEAST to obtain a mOTU-based tree for this analysis. We labeled each sample with their respective geographic distributions: A (Sumatra-North, comprises all samples from the Malacca Strait and Siam river systems including one from their neighboring watershed in the western slopes of Bukit Barisan), B (Sumatra-Central, comprises all samples from the North Sunda river system including one from their neighboring watershed in the western slopes of Bukit Barisan), C (Sumatra-South, comprises all samples from the East Sunda river system including one from their neighboring watershed in the western slopes of Bukit Barisan), D (Java), E (Borneo), F (Asia). We ran the analysis using the BayAreaLIKE + J model, which was identified as optimal during model selection by the program, set the number of maximum areas at each node to six, while leaving all other settings at default values.

Data availability
Data generated for this study is available in Supplementary information. Sequences are deposited in GenBank with accession numbers provided in Table S1.