Harrison’s rule corroborated for the body size of cleptoparasitic cuckoo bees (Hymenoptera: Apidae: Nomadinae) and their hosts

Harrison’s rule, that body size is positively correlated between parasites and hosts, has been reported in a range of taxa, but whether the rule is applicable to cleptoparasitic insects is poorly understood. Subfamily Nomadinae, the largest group of cleptoparasitic bees, usurp the nests of a variety of host bees. Within the subfamily, Nomada exploits the most diverse hosts, using at least ten genera from five families. Here, we reassess the phylogeny of Nomadinae, including the expanded sampling of the genus Nomada, to explore host shift fluctuations throughout their evolutionary history and test the applicability of Harrison’s rule for the subfamily. Our phylogenetic results are mostly congruent with previous investigations, but we infer the tribe Hexepeolini as a sister taxon to the tribe Nomadini. Additionally, the results reveal discrepancies with the traditional classifications of Nomada. Ancestral state reconstruction of host use indicates that, early in their evolution, parasites used closer relatives, before attacking less related groups later. Lastly, we confirm Harrison’s rule in Nomadinae, supporting that body size dynamics influence the host shifts of cleptoparasitic bees.

Body size influences many aspects of organismal biology, such as patterns of resource use 1-3 , habitat colonization potential 4,5 , and ecological strategy 6,7 . Among species with strong ecological interactions, body size can play an even stronger role, and this is especially important for parasites and their hosts 8 . The importance of size in parasitism is highlighted in Harrison's rule (HR), which states that host and parasite body sizes generally positively covary 9 .
Given that parasite body size depends on that of their hosts to some degree, we must explore the evolutionary trends of host-parasite association to test for morphological similarity and divergence between parasites and hosts. Although HR is a common pattern across the animal kingdom, it remains poorly understood from a macroevolutionary viewpoint; the rule has been mostly demonstrated using phylogenetically-independent comparisons to date (i.e. simple allometry between hosts and parasites 10-15 ,). Exploring HR under a molecular phylogenetic framework may offer new insights into the underlying adaptive basis and historical context 16 , while also providing invaluable quantitative metrics of relatedness. For example, are body size correlations a result of co-speciation or does it largely result from shifts to a host that has a similar body size? For such a study, it is ideal to focus on a monophyletic lineage of obligatory parasites with both high body size variation and good host-parasite association data.
Also, HR has been documented in many parasitic organisms: parasitic nematodes and their hosts 12 , parasitic barnacles and decapods 17 , avian lice and birds 18 , and fleas and rodents 19 . HR has been tested for several insect groups previously, in terms of both traditional parasitism and herbivory, but most insect groups remain unexplored 11 . One parasitic organism for which HR has never been extensively and quantitatively tested is the cleptoparasitic bees. About 13% of all bees and 20% of bees in the family Apidae exhibit the cleptoparasitic

Results
Phylogenetic analysis. Regarding taxon sampling, we targeted species that have both host and body size information without bias. In total, 106 species including 2 species for the outgroup and 104 species for the ingroup were used. We extracted novel sequences from 35 species, and 71 species were obtained from NCBI (Table S1). About 30% of ingroup taxa belong to Nomada, which has an relatively unstable phylogenetic position and possesses the most host information in the subfamily. The dataset used for phylogenetic reconstruction contained a total of 4590 bp (657 bp of COI, 709 bp of Ef1α, 1463 bp of Nak, 459 bp of Opsin, 843 bp of PolII, 459 bp of Wng). Phylogenies obtained through Bayesian inference (BI) from MrBayes and maximum likelihood (ML) from IQtree indicate strong support for the monophyly of Nomadinae (Bayesian posterior probability, PP = 100, Maximum likelihood bootstrap values, BS = 99) (Fig. 1).
Melectini was recovered sister to the remaining cleptoparasitic lineages with high support (PP = 100, BS = 98) (Fig. 1) The tribe Nomadini was also strongly supported as monophyletic (PP = 100, BS = 100), and it was recovered sister to Hexepeolini (PP = 78, BS = 94) (Fig. 1). In comparison to the tribal relationships, the infra-tribal relationships of Nomadini were more complex. Some simple cases where single species represented species-groups could be valid, for example the roberjeotiana species-group (N. roberjeotiana) and basalis species-group (N. japonica) (sensu Alexander & Schwarz, 1994 41 ). Similarly, three members of the furva group, N. distinguenda, N. discedens, N. okubira, were supported as monophyletic with strong support (PP = 100, BS = 100). However, the ruficornis species-group was more problematic. For example, N. lathburiana, which was previously designated as the ruficornis species-group 41 , was part of the bifasciata species-group. Also, N. aswensis which is supposed to be part of the ruficornis species-group was part of a clade with N. nipponica (tripsona species-group), N. armata, and N. ginran (armata species-group). A difference between BI and ML was found in the placement of N. alboguttata. In BI, the species formed subclade with N. fulvicornis jezoensis, N. fervens, and N. marshamella. On the other hand, it grouped only with N. fulvicornis jezoensis according to ML.
Ancestral state reconstruction of host use. The ancestral states of host use in Nomadinae were analyzed at the family level to enable easier interpretation of the results, avoid an overestimation of shifts, and accommodate program analytical limits. Except for a few nodes, parsimony and RJ-MCMC analyses were largely congruent. According to the Bayestraits analysis, the common ancestral host of subfamily Nomadinae was the family Apidae, with a probability of ~ 63% (Fig. 2). This host family is also inferred as the ancestral host for multiple lineages in Nomadinae. Additionally, we found that the host use of Apidae by Melectini, Ericrocidini, Ammobatoidini was conserved. On the other hand, frequent reversal host shifts between Apidae and Andrenidae were detected. After the first reversal host transition from Andrenidae to Apidae happened in the tribal combination of Isepeolini + Parepeolini + Protepeolini + Rathymini + Ericrocidini, host switching also occurred from Apidae to Colletidae and Melittidae (Fig. 2). The common ancestral host of Epeolini was reconstructed as Family Apidae with a probability of ~ 62%. However, a host switch to the family Colletidae in the genus Epeolus was observed with a probability of > 99%.
In the clade of Ammobatoidini + Neolarrini + Hexepeolini + Nomadini, their common ancestral host use was yielded as Family Andrenidae with a probability of > 96% (Fig. 2). Thereafter, a host shift occurred from Andrenidae to Halictidae at least three times in this group. To be specific, the first transition happened in the tribe Biasitini, and the remaining two transitions in different clades of the tribe Nomadini. In the tribe Nomadini, there was one switch from Andrenidae to Apidae within our sampling.
The most frequent host switches were observed in the tribe Ammobatini (Fig. 2), and the tribe showed a tendency of reversal host shifts from Andrenidae to Apidae. However, the ancestral host family reconstruction for the tribe was ambiguous and recovered associations with multiple families.
Allometry and ancestral state reconstruction of body size. The body length of the Nomadinae was on average 9.08 ± 0.36 (mean ± S.E.), widely ranging from 2.5 to 23 mm, and the intertegular distance (ITD) was 2.06 ± 0.11, ranging from 0.48 to 6.82 mm. The mean body length of the hosts was 10.66 ± 0.38, ranging from 3.6 to 23.28 mm and host ITD was 2.84 ± 0.17, ranging from 0.73 to 8.15 mm (see Supplementary Data 1). Our linear regression analysis strongly suggests that the size between Nomadinae and its hosts is highly related in terms of both body length (R 2 = 0.6879, P < 0.05) and ITD (R 2 = 0.7620, P < 0.05).
The ancestral state of ITD and body length was of moderate length (mean of full body length of Melectini: 11.06 mm, median ITD of Melectini: 3 mm; mean full body length of hosts: 13.43 mm; mean ITD of Melectini: 3.59 mm), and the length and width of body forms have evolved from medium to extremes (Fig. 3). It was observed in the tribe Ericrocidini that the body evolved to become longer and wider, but this phenomenon was found only in the tribe Ericrocidini. On the other hand, becoming shorter and narrower was observed in multiple lineages. This distinctly recognizable pattern of shrinking was detected in the tribes Ammobatini, Brachynomadini, Ammobatoidini, Neolarrini, and Nomadini.
Intriguingly, we found potential indicators that changing size is related to host shifting across the phylogeny. To be specific, a considerable body size increase was detected in the tribe Ericrocidini as the host changes from Andrenidae to Apidae. In contrast, both body length and width remained relatively similar within the genus Epeolus, which is parasitic on just the genus Colletes. This was not always the case, however, as seen among the species-groups of Nomada. For example, N. flavoguttata, which belongs to the Nomada ruficornis species-group and parasitizes the genus Andrena, showed a similar body size with species of the Nomada furva species-group even though that species-group specializes on a different genus, Lasioglossum.

Discussion
Most of the relationships recovered in this study correspond with a prior investigation that used phylogenomic sampling 31 31 . Third, the tribe Hexepeolini was inferred as the sister group of the tribe Nomadini. The placement of this tribe has varied even among molecular studies as sometimes it is thought to be sister to the Ammobatoidini, Ammobatini, or Neolarrini 26,28-31 . However, our results are relatively similar to those of the most recent study 31 , where Nomadini was sister to Hexepeolini + Neolarrini. The discrepancy between our results and previous Sanger studies is likely due to differences in taxon sampling, since most prior studies have used relatively few samples. In general, however, we find the most recent Nomadinae phylogenomic study is likely representative 31 , given the depth of data used and their good coverage across taxonomic groups.  www.nature.com/scientificreports/ Odanaka and colleagues recently conducted the first densely-sampled molecular phylogenetic investigation of Nomada using UCE phylogenomic data with 16 species-groups, and their investigation suggested that the odontophora species-group is sister to the rest 45 , contrasting with our results that suggested the roberjeotiana species-group was sister to the rest. Notably, however, we did not include a sample of the odontophora speciesgroup and consequently, here and elsewhere, it seems likely that the differences are due to taxon sampling. Further, difficulties in assigning species-group membership correctly are well known for this genus 40 , such that this could explain some disagreements between our two studies. For example, they recovered the superba and ruficornis species-groups as paraphyletic, while our study found that the armata species-group is paraphyletic, and ruficornis species-group is polyphyletic 45 . Our multi-gene phylogeny supports the designation of N. ginran as a member of the armata species-group as suggested by 42 . On the other hand, in disagreement with Alexander and Schwarz 41 , N. aswensis and N. kaguya (previously treated in the ruficornis species-group) and N. nipponica (previously tripsona species-group) were also recovered in the armata group. Additionally, N. lathburiana (of the ruficornis species-group) was recovered in the same clade with N. succincta, N. comparata, N. goodeniana, which belong to the members of the bifasciata species-group. Therefore, further investigation with increased taxon sampling is necessary to resolve the true composition and relationships of the Nomada species-groups. Again, differences in results between our study and the most recent investigation may be due to differences in taxon sampling, since we primarily focused on species with valid host information. In addition, Odanaka et al. examined Nomada phylogeny at a deeper level by using the UCE phylogenomic data with a much denser sampling pool 45 . Given that the phylogenomic approach is possible to provide dept and resolution to the data, the incongruencies between the latest study and the present work are presumably due to the different research methodologies.
Our results largely agree with prior studies, with the family Apidae as an original host and multiple switches to other groups and sometimes back to Apidae. In this way, nomadines appear to initially follow Emery's rule, with parasites only attacking close relatives as seen with Melectini attacking Anthophorinae, but later with many groups attacking entirely different families very successfully 31 . www.nature.com/scientificreports/ In comparision to the most recent reconstruction of Nomadinae 31 , our study determined that host shifting from Apidae to Andrenidae in Nomadinae occurs much earlier. This is because the tribe Caenoprosopidini is shown as sister to the Isepeolini + Epeoloidini + Protepeolini + Rhathymini + Ericrocidini subclade which belongs to the Ericrocidini line (sensu Sless et al., 2022) in our investigation. Conversely, in the previous study, Caenoprosopidini was placed within the Nomadinae line 31 . Denser sampling of these groups will help resolve these and other disagreements.
Host specificity to certain groups has been frequently used to define taxonomic units, including in cuckoo wasps, fig wasps, and some mites 24,[49][50][51] . Some similar examples exist in Nomadinae. For example, the host genus of Epeolus is well known to be limited to the genus Colletes 52 , and Ericrocidini is limited to Centridini 39 . The common feature of these two examples is that their hosts have specialized nesting strategies, and this may explain some instances of conservatism. Most colletid females apply a cellophane-like layer with their short, bilobed glossa in brood cells 39,53 . Conversely, many hosts of the tribe Ericrocidini collect floral oil for their brood-cell construction 20 . Applying floral oil and other secretions may make it difficult for parasites to access these nests 54 . Both of these specialized nesting behaviors, especially oils which may be hard to digest, could physically inhibit the offspring of parasites 55 . It may be that the adaptations for using these nests, in turn, make it harder for them to switch to other host groups with less specialized behaviors, further enforcing patterns of conservatism 56 .
Of course, exceptions to conservatism of host use exist, such as the use of many different hosts by Nomada. Coincidently, this group also represents another inconsistency with the most recent study 31 . We found three host shifts in the tribe, while they found one. There are over 800 species of Nomada 35 and many of them have unknown hosts, so this result depends heavily on taxon sampling. With greater sampling, we would expect additional host shifts. For instance, the host of the Nomada emarginata Morawitz (ruficornis species-group) is recorded as Melitta haemorrhoidalis Fabricus 57 , which belongs to the family Melittidae, but it was not included in our analysis. Similarly, the unused Nomada articulata Smith (erigeronis species-group) takes Agapostemon sericeus Foster and A. virescens Fabricius as their hosts 58 , and this could add another shift. As there are numerous species with unknown or unconfirmed hosts, many of which might invoke additional host shifts when added to the phylogeny, it is difficult to estimate the expected total number of host shifts for this group at this time.
Our size allometry analysis strongly supports the positive correlation between parasites and hosts, commonly referred to as Harrison's rule (HR 9 ;). Further, we found multiple instances where host shifts were linked with body size changes. It is likely that denser taxon sampling will reveal more such linkages between shifts in host use and body size, and it may be that this is a broader co-evolutionary pattern in Nomadinae.
Similar to prior studies on other taxa 59-62 , the body size of cleptoparasitic bees was typically smaller than their hosts in this study. We postulate that there is a "lock-and-key" relationship between the Nomadinae and their hosts (referred to in lice prior, based on the need to fit well on hair 13 ;). Most of the hosts of Nomadinae are soil-nesting, solitary bees. Females make their nests underground by digging, and their nest tunnels are typically just large enough to allow them passage. Because cleptoparasitic bees must enter nests to lay their eggs, they are constrained by the size of these tunnels, nest cells, and even the amount of food provisioned for offspring 25,63,64 . Our investigation supports this idea, with host shifts or stability seemingly influencing body size. For example, as shown in Fig. 3, there is relatively little size fluctuation in the genus Epeolus, all of which attack Colletes 52 . On the other hand, Nomada, which parasitizes a vast array of groups, shows scattered size-shifting corresponding with host switches. The ability of cleptoparasitic bees to modulate their size in response to hosts may contribute to the ability of these bees to switch hosts more readily, which could in turn enable adaptive radiations generated from switching to many new hosts (a possible component of the success of the huge group Nomada). Of course, plasticity in host recognition via visual and chemical pathways, as well as adaptability to different nectar and pollen resources, would also be important components.
Another example that highlights the relationship between host switching and body size modulation is found in the tribe Ammobatini. Even though the genus Melanempis was excluded from our study because its host is unknown, it is one of the largest nomadines and warrants discussion. According to the most recent study 31 , Ammobatini includes the small-sized Oreopasites (as small as 2.2 mm) and Melanempis (up to 22 mm), demonstrating huge size variation in the group 36,65 . Alongside this wide range of body sizes, they also parasitize many families: Apidae, Andrenidae, Colletidae, and Halictidae (See supplementary data in 31 ). With further data from this group on host associations and body size, and greater taxon sampling, it could be an ideal finer-scale study for looking at the influence of host switching on body size.
Body size is thought to be one of the factors that could contribute to cophylogenetic patterns between hosts and parasites 66 . Food resource requirements being similar among close relatives can help explain such patterns, as they should be relatively similar 13 . Indeed, the ancestral host of Nomadinae was revealed as Apidae in this study (specifically Anthophorinae 31 ,), so co-evolution may have taken place more closely early in the evolution of the nomadines. However, behavioral flexibility can empower parasites to invade novel hosts 56,67 . Such potentials are enhanced when alternative hosts have similar-sized nests 25 , or when their increased abundance makes it more likely for them to be encountered, especially if in the same microhabitat 68 . Host chemical cues also likely play an important role but are relatively poorly understood 69 . Such factors may help explain how single species can exploit many varied hosts 70 . For example, Nomada flavoguttata parasites various Andrena species such as Andrena falsifica, A. minutula, A. minutuloides, A. semilaevis, and A. subopaca 35 . Notably, all these host species belong to the same subgenus, Micrandrena, which is well-known to have diminutive body size, and may also have similar chemical cues and metabolic requirements 71 . However, occasionally, a single Nomada species exploits multiple genera or even families as hosts, such as in N. imbricata, parasitizing Andrena and Halictus 70 . With greater emphasis on natural history studies to make more host-parasite associations, even more disparate host uses might be discovered.
With the combination of our multi-locus phylogeny and allometric data, our understanding of the evolutionary relationship of size fluctuation between cleptoparasitic bees and their hosts has greatly improved. Based on www.nature.com/scientificreports/ the scenarios and examples discussed above, host switches appear closely related to changes in body size. The ability to adapt to new hosts more quickly in this way could increase speciation rate in the presence of potential new hosts, explaining highly diverse groups like Nomada. Ultimately, we confirm Harrison's rule, the positive correlation of body size between cleptoparasitic bees and hosts, and with further sampling we expect to find an even stronger impact of body size on the evolution of this group.

Methods
Taxon sampling. In total, we used 104 species for the ingroup. For the outgroup, Apis cerana and Colletes compactus were selected to represent a range of relatedness and performed well in pilot analyses, compared to the relationships recovered in a recent phylogenomic analysis 31 . Sequences for 71 species from published papers were taken from NCBI, and we newly added 35 species that have not been analyzed in previous phylogenetic investigations (Table S1). To minimize missing data and ensure better resolution in the phylogenetic tree, we used sequence data only if the number of available genes was greater than or equal to three. The aim of the sampling was to use species that have both host and body size information without bias. When species did not fulfill this condition, they were excluded from this study. However, we made exceptions where necessary to ensure that the resultant phylogenetic reconstructions were reliable, with reference to prior works informing our choices 31,72 .
We placed special focus on sampling the genus Nomada as it has highly-diverse host records 31 . Specifically, we expanded the sampling of the furva species-group because they are distinctly small-bodied and the parasite of genus Lasioglossum (Hymenoptera: Halictidae). Two species, N. aswensis, and N. kaguya, which belong to the ruficornis species-group but show similar biology with the furva species-group, were also included in this study. We included two type species, N. ruficornis from ruficornis species-group and Nomada armata from armata species-group, to fully investigate the composition of these species-group concepts by Alexander 40 , although greater sampling is necessary to fully define the species-groups 45 . The study used DNA samples of Finnish Nomada species from the Genomic Resources collection (https:// laji. fi/ en/ theme/ luomu sgrc/ instr uctio ns) of the Finnish Museum of Natural History Luomus (https:// www. luomus. fi/ en), and the HTTP-URIs of the DNA samples are listed in Table S1. In total, our sampling of the genus Nomada represents seven of the 16 speciesgroups, including a range of host information in three families representing Apidae, Andrenidae, and Halictidae. We could not represent all family-level hosts for this group because there are few species attacking some groups, such as Melittidae, and they were more difficult to sample. DNA extraction, PCR amplification, and sequencing. Total genomic DNA was extracted by grinding up detached midlegs or heads of either alcohol vouchers or dried specimens; forceps used for extraction were sanitized before by flame and rinsed in 99% EtOH between each specimen. DNeasy Blood and Tissue kit (QIA-GEN, Inc.) was used for the DNA extraction according to the manufacturer's protocols and stored in − 20 °C. DNA vouchers were deposited in the Insect Biosystematics Laboratory, Seoul National University.
PCR products were amplified using Accupowder PCR Premix (Bioneer, Daejeon, Korea) in a 20 µl reaction mixture. In case DNA extraction was needed from bad quality specimens, we redesigned the partial primers in threefold. The primers used in this study and specific PCR conditions are listed in Table S2. After the amplification process, PCR products were purified and sequenced at Bionics CO. (Seoul, Korea).
Using the SeqMan Pro version 7.1.0., raw sequence data were assembled, checked, and trimmed. Sequence alignment of six genes were performed in MAFFT version 7 (https:// mafft. cbrc. jp/ align ment/ server/) 79 . All the sequences were adjusted in Mega 7 with the amino acid translation option. The aligned sequences were combined using SequenceMatrix Windows ver. 1.8 80 .

Phylogenetic analyses. Phylogenetic analyses were performed using both Bayesian inference (BI) and
Maximum likelihood (ML). The Bayesian inference analysis was conducted with MrBayes 3.2.7a 81 . The best substitution model for BI was selected for each partition under the Bayesian Information Criterion (BIC) using ModelFinder in IQTREE and the same protocol was used in searching for the best substitution model for ML 82,83 . GTR + F + I + G4 for COI, TIMe + I + G4 for ef1α, TIM3 + F + I + G4 for Nak, TIM3e + I + G4 for Opsin, TIM2 + F + I + G4 for PolII, and TIM + F + I + G4 for Wng were selected as the best fitting models for BI. For ML, GTR + F + I + G4 for COI, TIM + F + I + G4 for ef1α, TIM3 + F + I + G4 for Nak, TIM3e + I + G4 for Opsin, TIM2 + F + I + G4 for PolII, TIM + F + I + G4 for Wng were selected 84 . However, since TIM models cannot be used in MrBayes, GTR + I + G was used to run the BI analysis except COI. For the MrBayes analysis, we ran 20 million Markov chain Monte Carlo (MCMC) generations and trees were sampled every 100 generations. We ran one cold chain and three heated chains for each MCMC analysis. The first two million sampled trees were discarded as burn-in. Branch support for the maximum likelihood analysis was computed using the UltraFast bootstrap approximation (UBS 85 ;) with 1,000 replicates.
Size allometry. We measured two traits of the Nomadinae and their hosts. One is the intertegular distance (minimum distance between the tegulae), which is a useful estimator for the size of bees 86 . The other is total body length, which was measured as the maximum length distance from the head to end of the final tergite exclusive of exserted stingers or genitalia. (if curled, this was accounted for using a multi-stage measurement). In total, we obtained about 1,300 body size data from Nomadinae and their hosts (Body length: 73; ITD: 65 out of 106 species of Nomadinae and their corresponding hosts). Among them, we used specimens of 24 species (host: 7, www.nature.com/scientificreports/ Nomadinae: 17) and prepared photographs using the Microscope (DM 4000B, Leica Microsystem, Wetzlar, Germany) with a USB digital camera (Infinity3, Lumenera Corporation, Ottawa, Ontario) and measured the value of the traits using the measurement option. Since there can be a size difference between male and female bees, measurements were only taken from females. We also collected data from literature such as taxonomic papers, online-accessible specimen photographs in museums, or from taxonomists (see Supplementary Data 1). The number of measured individuals was different per species because we gathered the data from various sources. Therefore, we used the average value of measurements for the data analysis. Linear regressions were conducted using SPSS Statistics 25 (IBM, Armonk, NY, USA) to gauge the dependence of nomadine body length on host body length, as well as nomadine ITD on host ITD.
Ancestral character state reconstruction. The host information of each nomadine species was extracted from the literature (See supplementary data 1). Where multiple species were reported as hosts, the data was prioritized when the parasitic larvae were found in their nests or when direct intrusion of the parasite to the host nests was observed. However, when such confirmation was absent, data arose from phenological synchrony between parasites and hosts and information provided by taxonomists, including hypotheses raised in the literature. The hosts of Nomadinae were coded as five discrete states by families: (A) Andrenidae, (B) Apidae, (C) Colletidae, (D) Halticidae, and (E) Melittidae. This may underestimate the number of shifts because different lineages are used in some families, but as it is difficult to exactly know when switches occurred this is in some cases conservative, and avoids the artificial inflation of the number of shifts. The trace character history function in in Mesquite 3.31 was used to map the evolutionary history of host use on a single tree with parsimony approach 87 was adopted. The probability of the ancestral state of each node was calculated by Bayestraits v3.0 88 using reversible jump Markov Chain Monte Carlo (RJ-MCMC). An exponential distribution was implemented, seeding from a uniform prior in an interval of 0-100. We ran 50 million iterations, sampling every 1,000th iteration. The first million iterations were discarded as burn-in. Acceptance rates were automatically adjusted and achieved in the preferred range of near 35%.
To trace the evolutionary history of intertegular distance and body size, we used ancestral state reconstruction using the maximum likelihood (ML) method in the R package 'Phytools' 89 . The resulting Bayesian tree was converted into a dendrogram, the missing data vector was estimated using ML, and the ContMap function was used to visualize their evolutionary history.