Reticulate evolution in the Pteris fauriei group (Pteridaceae)

The Pteris fauriei group (Pteridaceae) has a wide distribution in Eastern Asia and includes 18 species with similar but varied morphology. We collected more than 300 specimens of the P. fauriei group and determined ploidy by flow cytometry and inferred phylogenies by molecular analyses of chloroplast and nuclear DNA markers. Our results reveal a complicated reticulate evolution, consisting of seven parental taxa and 58 hybrids. The large number of hybrid taxa have added significant morphological complexity to the group leading to difficult taxonomic issues. The hybrids generally had broader ranges and more populations than their parental taxa. Genetic combination of different pairs of parental species created divergent phenotypes of hybrids, exhibited by both morphological characteristics and ecological fidelities. Niche novelty could facilitate hybrid speciation. Apogamy is common in this group and potentially contributes to the sustainability of the whole group. We propose that frequent hybridizations among members of the P. fauriei group generate and maintain genetic diversity, via novel genetic combinations, niche differentiation, and apogamy.

www.nature.com/scientificreports/ hybridization. For example, Kuo 53 suggested that Pteris wulaiensis C. M. Kuo was a hybrid between P. faurei Hieron. and P. bella Tagawa. The Pteris fauriei group is mainly distributed in Eastern Asia and taxa have diverse geographic distributions and niche fidelities. Pteris fauriei, P. arisanensis Tagawa, P. biaurita L., and P. latipinna Y.S.Chao & W.L.Chiou are common and distributed widely 51,54 , but most taxa in the P. fauriei group are endemic to small regions. For example, P. boninensis H.Ohba, P. laurisilvicola Sa.Kurata, P. natiensis Tagawa, P. satsumana Sa.Kurata, and P. yakuinsularis Sa. Kurata are limited to Japan 55 , P. confusa T.G.Walker is found in India and Sri Lanka 56 , P. austrotaiwanensis Y.S.Chao is only reported in south Taiwan 50 , and P. kawabatae Sa.Kurata, P. minor (Hieron.) Y.S.Chao, and P. wulaiensis are recorded in Japan and Taiwan 51,57 . In regards to niche diversity, most taxa occur in shady environments, like P. latipinna, or semi-open habitats, like P. fauriei, although P. minor occurs in open habitats 51,58 . We addressed the following questions of the P. fauriei group: (1) Are any species of the P. fauriei group of hybrid origin? Chloroplast and nuclear DNA data were analyzed phylogenetically and integrated with ploidy level data to infer a possible network of hybridization. (2) Is morphological diversity related to hybridization? Morphological characteristics and genetic composition of each taxon were compared. (3) Is niche diversity related to hybridization? General attributes of ecological niches were compared among taxa. (4) What is the role of apogamy in the evolution of the P. fauriei group? Contributions and/or disadvantages of apogamy in hybridization were also discussed.

Results
Morphological and habitat characters. Samples of the Pteris fauriei group were collected from Eastern Asia, including Japan, Taiwan, China, Vietnam, and Malaysia, and included 20 taxa (18 species; Supplementary Table 1). No material of P. kiuschiuensis var. centrochinensis was collected for molecular analysis, and P. confusa and P. oshimensis var. paraemeiensis were amplified only for cpDNA markers. Some sampled plants morphologically matched the type specimens (Table 1); other plants were attributed to the most morphologically similar species. We found that the taxa had different niche attributes, such as light intensity and elevation ( Table 1). All species occurred in shady or semi-shady environments, except P. minor, which was found in open and sunny areas. Most taxa grew at low elevations (≤ 1000 m), but P. arisanensis, P. biaurita, and P. setulosocostulata were located at both low and higher elevations (> 1000 m).
Ploidy levels and reproductive modes. There were diploids and triploids in the P. fauriei group based on flow cytometry data; some named taxa had both diploids and triploids. Previous studies on ploidy levels in the P. fauriei group are cited in Table 1. There were seven taxa with homozygous nDNA genotypes that could be potential parental taxa (inferred by nDNA described in the section below): P. arisanensis, P. biaurita, P. boninensis, P. latipinna, P. minor, P. oshimensis var. oshimensis, and P. wulaiensis. These are all diploids except for P. arisanensis, which has some homozygous triploid individuals. Only three of the diploid taxa-P. boninensis, P. minor, and P. oshimensis var. oshimensis-had 64 spores per sporangium, which indicated that they reproduce sexually. The other plants had 32 spores per sporangium, which we inferred to be apogamous (Table 1; Supplementary Table 1). It has been reported that sexual Pteris species produce 64 spores per sporangium whereas apogamous species produce 32 spores per sporangium, especial in P. fauriei and P. minor 59 . Each plant produced only one type of sporangium.
The phylogeny of chloroplast DNA and nuclear DNA. The sequences of rbcL, matK, IBR3, and Knox3 of the P. fauriei group in this study (Supplementary Table 1) were clustered and named as haplotypes and allele types. Phylogenetic statistics are shown in Table 2. There were 20 cpDNA haplotypes; some were specific to one species and some were shared by several species. We included different taxa, different species and/or one species with different cpDNA haplotypes, of the P. fauriei group in the cpDNA phylogenetic analyses. For each species, samples with different haplotypes were included, and the samples with morphology identical to type specimens were marked (Fig. 1). The cpDNA topology was similar to the previously published Pteris phylogeny 60 , which resolved the P. fauriei group as part of sect. Campteria, although the supporting value was not very high.
The phylogenetic analysis of the Knox3 alleles supported six clades, which we labeled B, C, D, E, F + G, and H (Fig. 2a.). Group A is related to P. minor, and groups F and G correspond to P. arisanensis and P. biaurita, respectively, and are interdigitated. The IBR3 topology (Fig. 2b) approximately corresponded to that based on Knox3. Together with ploidy data, the genotypes of IBR3 and Knox3 were determined (Supplementary Table 1). Most of the genotypes were found in two or more samples from each taxon, except P. kawabatae, in which each individual sampled had a unique genotype ( Table 1). The IBR3 gene exhibited fewer alleles than the Knox3 gene (Table 1) (Table 1). We also phased the alleles of the two nDNA genes, to infer subgenome evolution of this group (Fig. 3). There is a mean posterior probability of − 3437.35, and the effective sample size (ESS) was 292. The results support that the two nDNA genes present similar evolutionary histories, and the inferred phylogeny is generally consistent with the pattern of two DNA phylogenies in Fig. 2. Potential parental taxa with homozygous nDNA genotypes were shown empty in the nDNA columns. Most of the posterior probabilities of the phase assignments were high (heatmap in blue), and potential missing data (not be amplified alleles) are shown as having low posterior probability (heatmap in red).
The following analyses were based on the Knox3 gene results mainly because of their higher resolution than the IBR3 gene and show clearer lineages than the phasing patterns to infer hybridization relationships between potential paternal taxa and hybrids. Subsequently, a phylogenetic network was inferred based on Knox3 alleles and genotypes (Fig. 4). The crossing patterns (taxa in Table 1) were inferred. Lineages of A to G groups were colored separately, and then the lineages of putative hybrid taxa were indicated. The topology presents an evolutionary radiation involving hybridization; many hybrid taxa arose from only a few parental lineages. Parent and hybrid taxa assignment. By comparing cpDNA haplotypes ( Fig. 1) and the phasing of nDNA genotypes (Fig. 3), the maternal lineages of the nDNA genotypes were inferred (Table 1). For example, the cf, cf ' , and cb cpDNA haplotypes were always present along with the group A allele of nDNA, and ca of cpDNA was found with allele D7 of Knox3 nDNA. A reticulogram based on the Knox3 gene was constructed, onto which we mapped ploidy levels, habitats, and reproductive modes (Fig. 5); the taxa correspond to the taxa in Table 1. Sexual plants were few, so we assumed that the probabilities of backcrossing with parents and introgression were low, especially with the apogamous diploid. Furthermore, the two nDNA markers presented similar topologies.
Based on the Knox3 gene, 72 genotypes were inferred, including 14 putative parental genotypes and 58 genotypes of hybrid origin (Table 1). All species were identified based on morphology. Plants with homozygous genotypes or some sexual diploid genotypes possessing alleles from the same groups, are proposed as possible parental taxa. Seven putative parental taxa were discovered. Strictly speaking, only samples morphologically identical with type specimens of P. arisanensis (F group, F7F7F7), P. boninensis (A7A7), P. minor (A group except A7 & A4, A1A1, A1A13, A1A6, A6A6), and P. latipinna (D7D7) were homozygous, which suggests that they were not hybrid taxa. Most taxa in the P. fauriei group, however, appeared to be of hybrid origin. Although the other three putative species had homozygous samples, such as P. biaurita belonging to group G (samples with G1G1), P. oshimensis var. oshimensis belonging to A4 (Kuo 3445, A4A4), and P. wulaiensis belonging to D4 (Ebihara et al. 3234, D4D4), they also had heterozygous taxa. It was difficult to separate the homozygous samples and heterozygous samples based on their morphology.  www.nature.com/scientificreports/ Pteris latipinna (D7D7) was the putative parental lineage of 18 hybrid taxa (Fig. 4), far more than any other species. Pteris wulaiensis (D4D4) was the second most important parental lineage, contributing to eight taxa. Homozygous P. oshimensis var. oshimensis (A4A4) was inferred to be the maternal parent of hybrid P. oshimensis var. oshimensis (A4D1) and P. cf. fauriei (A4D7).
Genotypes resolving puzzles of morphology and habitats. The diverse morphology of the P. fauriei group appeared to be an outcome of a large number of hybridizations. By comparing morphological characters 51 and molecular data (Table 1), the association of specific morphological character states with particular genotypes was explored. Some character states were only found in individuals with particular alleles (Table 3). We present apparent associations between morphological and ecological variation and genotypic markers as hypotheses that could be tested with more detailed analyses. We found that hybrid taxa exhibited the morphology associated with their parents (Table 1).
Plants with homozygous genotypes provided more apparent evidence of character state/genotype associations, than did heterogeneous taxa (Tables 1 and 3). Number of pairs of basiscopic secondary pinnae were useful key characters in this group. In P. minor individuals with two or more pairs of basiscopic secondary pinnae or with tripartite laminae possessed alleles A1, A6, and A13. Pteris latipinna (D7D7) had the largest pinnae and the fewest lateral pinnae of all taxa; other taxa with wide pinnae had allele D7, including P. fauriei, P. kawabatae, P. laurisilvicola, and P. natiensis. It appears that the P. latipinna genome is associated with broad pinnae. www.nature.com/scientificreports/ The triangular costal veins of P. arisanensis are more-or-less intermediate between areolate costal veins, such as in P. biaurita, and the free venation of some other species, suggesting a possible hybrid origin of P. arisanensis. Homozygous P. arisanensis (F7F7F7, Chao2135), however, had triangular costal veins, and homozygous P. biaurita (G1G1, Chao3056) had areolate costal veins, thus hybridization does not appear to account for the areolate veins of the former species. By contrast, hybrids (F7G1*, F16G1) between two species exhibit areolate costal veins and were identified as P. biaurita.
It is thought that putative hybrids had intermediate phenotypes between two parental taxa but others were more similar to one or the other parent 67 . Hybrids with different genotypes could be morphologically similar if they had the same allele, from the same parental taxon. For example, P. kawabatae, P. kiuschiuensis, P. natiensis, and P. satsumana have basal pinna segments connecting to midribs and they share allele C3. Although no sample had a homozygous genotype of C3, the C3 allele was associated with this character. The segments' angle against midribs in P. kiuschiuensis var. kiuschiuensis and P. setulosocostulata was 80°-85°, wider than all the other taxa, and both taxa had H type alleles. Thus the H allele group appeared to be associated with the larger angle of the segments against midribs.
Most plants in the P. fauriei group are found below 1000 m elevation, except some plants of P. arisanensis (D7F4F5, F7F7F7) and P. biaurita (F7G1*, F16G1) can occur up to 2000 m. The Knox3 alleles of the F and G groups appeared to be associated with elevational range. Homozygous samples (i.e., putative parents) of P. biaurita (G1G1) were located in Southeast Asia at low elevation, whereas F group alleles were associated with occurrence at high elevation, above 1000 m. Hybrids (with both F and G group alleles) had broad elevational ranges (P. biaurita F7G1*, F16G1, low & high; Table 1, Supplementary Table 1).

Discussion
In recent evolutionary radiations, genetic and morphological divergences are low because of the short evolutionary time scales involved [68][69][70] . According to Chao, et al. 60 , the main fauriei clade arose around 5 Myr, and most taxa diversified around 2 Myr. Our cpDNA tree topology is based on the same markers as that study and supported a recent radiation of the P. fauriei group because most taxa were separated by distinctly shorter branches, corresponding to short time periods inferred in the study of Chao, et al. 60 . We revealed a complex hybridization network in the P. fauriei group, composed of morphologically similar taxa but also putative hybrid taxa exhibiting morphological intermediacy between probable parental species. This large number of hybrid taxa adds significant morphological complexity to the group leading to difficult taxonomic issues 51 . Furthermore, the evolution of the P. fauriei group accounts for a significant amount of the diversity of sect. Campteria in Pteris, i.e., 18 species, and especially the apogamous elements in Asia 45 . Sect. Campteria consisting of more than 60 species, could be the biggest section in Pteris 71 and has more hybridizations, such as in the P. quadriaurita complex 41 . Hybridization is a key mechanism facilitating the diversification of sect. Campteria.
Among the 18 species of the P. fauriei group, some species are phylogenetically distinct from the main clade. We have clarified the hybridization patterns of P. arisanensis and/or P. biaurita. Although P. confusa was thought to be a synonym of P. arisanensis 56 , our data suggest that they are different species 51 . We found that P. confusa and P. praetermissa have the same material lineage. However, because of limited materials, the hybridization involving P. confusa 7 is still unclear.
The assignment of homoeologs is complicated in taxa involved in hybridization and polyploidy [72][73][74] . Recently, more applicable methods for phasing gene copies into polyploid subgenomes are now in development, focusing on several loci 75 or target capture data 76 . Even though we found no evidence of possible introgression in the P. fauriei group, we revealed that multiple hybridizations contributed to offspring exhibiting diverse phenotypes. Frequent hybridization with low introgression is found in hybrid zones of cacti, apparently due to postzygotic isolation 77 .
Diversification of the P. fauriei group appears to be related to different environments. For example, light intensity of P. fauriei (apogamous; semi-shady), is an intermediate habitat between its parental taxa P. minor (sexual; open) and P. latipinna (apogamous; shady). By contrast, homozygous plants of P. arisanensis and P. biaurita were located at high and low elevation, respectively, and putative hybrids of the two species showed a wider elevational range than the parents, essentially representing the addition of ranges of both parental taxa.
Cases of apparently extinct diploid parental taxa are common in many fern complexes 9,48,49,78 and in the P. fauriei group as well. The distributions of extant parental taxa are much narrower than those of the hybrids, which might imply the decline of parental taxa. For example, P. latipinna is an important parent, while the known populations are fewer than ten in China and Taiwan and P. minor is limited to Iriomote Isl. (Japan) and Taiwan 51 . Their hybrid P. fauriei occurs across a broader area from Japan, Taiwan, and eastern China. Pteris oshimensis var. oshimensis is distributed across Japan, but the paternal plants are only found in Amami Is, Japan 51,52 . Thus, the ranges of hybrids often exceed those of their parents. The wide distribution of the P. fauriei group also suggests that hybridization has enhanced range expansion of the entire group 79 .
After long-distance dispersal, apogamy and gametophytic selfing could be adaptive because either is more likely to produce a new population than would sexual reproduction between gametophytes 80 . Apomictic angiosperms (i.e., asexually reproducing via seeds) tend to predominate in environments unfavorable for sexual reproduction, such as at higher latitudes and elevations, and have a wider distribution than their sexual relatives 81,82 . The apogamous triploid fern Myriopteris gracilis has a wide distribution which might be related to its ability to reproduce asexually 83 . Niche differentiation of apogamous and sexual fern taxa has also been documented previously 38 www.nature.com/scientificreports/  Table 1. ◂ those differences result from reproductive modes and/or higher fitness of certain genotypes needs to be explored further. Within the P. fauriei complex, apogamy is the main reproductive mode with only three sexual taxa. In general, asexual taxa, and lineages with low genetic variation, are thought to be less "adaptable" than are sexual taxa 42 . Furthermore, without recombination, deleterious mutations could accumulate and decrease the fixation of beneficial mutations (Muller's ratchet) 86,87 . However, new genetic combinations from hybridization could induce phenotypic differentiation of a hybrid from its progenitors 21,88,89 . In the P. fauriei group, the extant hybrid taxa, even being apogamous, are apparently the survivors of natural selection and have putatively retained adaptive genotypes that arose from hybridizations, even though some of their parental taxa may be extinct.
Furthermore, apogamous taxa have evolutionary advantages that might overcome the potential deleterious effects of low genetic variation. Epigenetic modification has been proposed to explain high fitness of apogamous taxa 34 . Moreover, Klekowski 90 proposed that genetic segregation is possible via homoeologous chromosome pairing during sporogenesis leading to genetic variation in apogamous ferns as was found in the study of Cyrtomium fortunei (Dryopteridaceae) 91,92 .
Of the seven diploid putative parental taxa, five are apogamous. We hypothesize that some parental taxa with apogamy in the P. fauriei group were primarily sexual and then apogamy developed subsequent to the production of hybrids. Apogamous homozygous diploids are more likely derived from sexual diploids rather than from hybrid taxa. The exact mechanisms of sexual taxa giving rise to apogamous taxa could be related to environmental factors such as light, water, and sugar [93][94][95] . Some critical genes expressed during reproduction could control the process of apogamy [96][97][98][99][100] . Apogamous gametophytes could contribute functional male gametes, which could account for the origin of some triploid taxa in the P. fauriei group.
The putative reticulate evolution of the P. fauriei group demonstrates a pattern of hybridization in evolutionary radiation, with limited or no introgression. Parental taxa have close phylogenetic relationships and novel genetic combinations in hybrids led to a rapid increase of phenotypic variation, such as diversification of niche fidelities, compared to parental taxa. Apogamy leads to the genetic fixation of some successful hybrid lineages under natural selection. While some parental species could be extinct, some of their genetic diversity will be preserved within the hybrids. Hybridization could lead to abundant taxa, more phenotypic variation, and broader distributions (as discussed above). Our ongoing and future research will explore the relationship of niche differentiation to hybridization, especially by exploring ecophysiological characters and environmental conditions. We hope to know if hybridization plays a role to maintain genetic variation and increase adaptation of ferns.

Methods
Sampling. Voucher specimens of samples were deposited at the herbarium of the Taiwan Forestry Research   Institute (TAIF; Supplementary Table 1). All sampled plants were identified based on the morphology of type specimens and original protologues; the plants that were morphologically the most similar to the types were marked as representatives of specific species. To identify the plants, morphological characteristics were examined following Chao, et al. 51 . Use of plant material in the study complied with relevant institutional, national, and international guidelines and legislation. Ploidy analysis and reproductive systems. Ploidy levels of samples were determined by flow cytometry, if fresh leaves were available, following the methods Chao, et al. 9 . Chromosome numbers of some samples had been counted in previous studies 46,52 . Previously published accounts of ploidy levels and reproductive modes of the P. fauriei group were also used ( Table 1).
The fact that the sexual Pteris plants produce 64 spores per sporangium and apogamous plants produce 32 or fewer spores per sporangium was used to infer the reproductive mode of each plant, corresponding to sexual or apogamous 29,31,59,101 . Five mature sporangia were picked randomly from each plant and spores were counted under a microscope.
DNA isolation, amplification. Materials for molecular analyses were preserved in silica gel. Total genomic DNA was extracted using a modified cetyltrimethylammonium bromide (CTAB) method 102 . Two plastid gene markers, rbcL and matK, were amplified by the primers employed in previous Pteris studies 57 . The two nuclear DNA (nDNA) putative single-copy genes, Knox3 and IBR3, were analyzed by single-strand conformation polymorphism (SSCP; see description below) and sequenced on the Illumina Miseq platform. Five primers for Knox3 were designed based on transcriptome sequencing data 103  www.nature.com/scientificreports/ than 500 bp. Both forward and reverse primers for Knox3 and IBR3 were synthesized with an 8-base barcode to produce amplicons; the barcodes were designed as recommended by Roche 105 .  Table 1. www.nature.com/scientificreports/ Single-strand conformation polymorphism. Some Knox3 data were analyzed by SSCP. The SSCP analyses were conducted by isolating the nDNA loci from PCR products for each individual through separation on SSCP gels, following the methods of Ebihara, et al. 8 . Individual bands were isolated from gels and purified by the Gel/PCR DNA Fragment Extraction Kit (Geneaid Biotech Ltd. Taipei, Taiwan) and amplified by primers Knox3PF2NE and Knox3PR3 for further Sanger sequencing.

Library preparation, sequencing and quality assessment for NGS. Each amplicon was extracted
with the illustra GFX PCR DNA and Gel Band Purification Kit (GE, UK). The amplicon mixtures were pooled with equal quantities of DNA, then a total amount of 150 ng of amplicons was used as input material for the DNA library preparations. The sequencing library was generated using the Truseq Nano DNA HT Sample Prep Kit (Illumina, USA) following the manufacturer's recommendations, and index codes were added to each sample. DNA fragments were ligated with the adapter for Illumina sequencing, followed by further PCR amplification.
Then the PCR products were purified (SPRIselect reagent, Beckman) and DNA size spectra were determined using an Agilent 2100 Bioanalyzer and quantified with a Qubit fluorometer (Invitrogen, Carlsbad, California, USA.). Finally, the DNA libraries were sequenced using the Illumina Miseq platform, and 300 bp paired-end reads were generated. Sequencing output was deposited in GenBank (Supplementary Table 1).
Bioinformatic analysis for nDNA data. Raw data were cleaned by FastQC 106 and Trimmomatic 107 to remove the adapter and low quality bases and reads. The forward and reverse reads were merged using PEAR 108 .
The cleaned data were demultiplexed and clustered (including chimeras removed), and chimeras were removed using PURC 109 . We kept the clusters comprising the three (for diploid samples) or four highest number of reads (for triploid or ploidy unknown samples). The barcode sequences were identified and removed, that is, each read was assigned to its source accession, then primers were trimmed. The sequences obtained from SSCP and Sanger  www.nature.com/scientificreports/ sequencing were used as the reference sequences for the clustering of consensus sequences. The sequences representing clusters of less than 100 original reads were removed. The consensus sequences of each sample were the alleles of nDNA genes and were ready for downstream phylogenetic analyses. The sequences of the Knox3 gene, both from SSCP (and Sanger sequencing) and NGS, were pooled together.
Data analysis. Sequences were automatically aligned using MUSCLE 110 and then manually edited with BioEdit 7.1.3 111 . Using DnaSp 112 , the sequences of the cpDNA and two nDNA markers were clustered by haplotypes and type of alleles, respectively. To exclude errors in amplicon data, we only kept the alleles found in two or more samples. We hypothesized that the taxa that had only one allele type were homozygous and classified them as parental taxa (non-hybrid taxa), regardless of whether their reproduction modes were sexual or apogamous. Allele numbers of genotypes were supported by ploidy data. When the allele number was smaller than the ploidy level, it was difficult to infer the exact allelic composition; in those cases, "*" stands for unassorted alleles 8,9 . Maximum likelihood analyses were performed for cpDNA haplotypes and nDNA sequences using the program GARLI v.2.0.1019 113 . Pteris species in sect. Campteris were included, and species in other sections, including P. bella Tagawa, P. formosana Baker, P. kidoi Sa.Kurata, P. longipes D.Don, and P. wallichiana J.Agardh, were chosen as outgroups 60,71 . For the three sequence datasets-cpDNA, Knox 3, and IBR3-three phylogenetic analyses were conducted by ten independent runs, from different random sequence-order starting trees, based on automatic termination following 10,000 generations without a significant topological change. The ML bootstrap support for each clade was assessed by performing 1000 bootstrap replicates, each replicate with one single tree search with the same search parameters as above. A 50% majority rule consensus tree was then calculated using PAUP* v. 4.0b10 114 . Gaps were treated as missing data. Then the nDNA genotypes, Knox3 and IBR3, were compared visually to the cpDNA markers to infer possible hybridization patterns. A phylogenetic network was visualized with Dendroscope 3 115 .
We also used homologizer 75 to phase the alleles of the two nDNA genes (linking the tree of the two genes) and infer a phylogeny. The potential parental taxa with homozygous nDNA genotypes were fixed firstly to infer the pattern of the subgenome (or allele) evolution of this group. The two genes were modelled by an independent Table 3. The genetic lineage inferences of morphological and habitat characters 51 , based on the Knox3 gene marker (Table 1) of Pteris fauriei group. The taxa of Pteris fauriei group exhibit the morphology corresponded to their own genotypes. The alleles corresponding to the character state are difficult to infer. www.nature.com/scientificreports/ GTR substitution model and exponential priors (mean = 0.1) on branch lengths. The analysis was based on a Bayesian Markov chain Monte Carlo (MCMC) approach for four runs; the run with highest posterior (mean) was selected. Trees were run for 15,000 generations and the first 0.1 trees were discarded as burn-in. We assessed convergence by Tracer v.1.7 116  www.nature.com/scientificreports/