Joëlle Barido-Sottani
email: joelle.barido-sottani@m4x.org

Alexander Pohle

Kenneth De Baets

Duncan Murdock

Rachel C M Warnock

Putting the F into FBD analysis: tree constraints or morphological data?

The fossilized birth-death (FBD) process provides an ideal model for inferring phylogenies from both extant and fossil taxa. Using this approach, fossils are directly integrated as part of the tree, leading to a statistically coherent prior on divergence times. Since fossils are typically not associated with molecular sequences, additional information is required to place fossils in the tree. We use simulations to evaluate two different approaches to handling fossil placement in FBD analyses: 1 using topological constraints, where the user specifies monophyletic clades based on established taxonomy, or using total-evidence analyses, which use a morphological data matrix in addition to the molecular alignment. We also explore how rate variation in fossil recovery or diversification rates impacts these approaches. We find that the extant topology is well recovered under all methods of fossil placement. Divergence times are similarly well recovered across all methods, with the exception of constraints which contain errors. We see similar patterns in datasets which include rate variation, however, relative errors in extant divergence times increase when more variation is included in the dataset, for all approaches using topological constraints, and particularly for constraints with errors. Finally, we show that trees recovered under the FBD model are more accurate than those estimated using non-time calibrated inference. Overall, we show that both fossil placement approaches are reliable even when including uncertainty. Our results underscore the importance of core taxonomic research, including morphological data collection and species descriptions, irrespective of the approach to handling phylogenetic uncertainty using the FBD process.

Introduction

Time-calibrated trees provide a crucial basis for hypothesis-testing in the life and earth sciences.

Phylogenetic dating combines molecular and fossil evidence, allowing us to reconstruct a timeline of events that are otherwise not directly observable. Within a Bayesian framework, temporal evidence is incorporated via the tree prior or tree model. The fossilised birth-death (FBD) process explicitly combines the lineage diversification and fossil recovery processes, providing an ideal model for inferring phylogenies from both extant species and fossil specimens [START_REF] Stadler | Sampling-through-time in birth-death trees[END_REF]. Using this approach, fossils (with or without character data) are directly considered as part of the tree [START_REF] Heath | The fossilized birth-death process for coherent calibration of divergence-time estimates[END_REF][START_REF] Gavryushkina | Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration[END_REF][START_REF] Gavryushkina | Bayesian total-evidence dating reveals the recent crown radiation of penguins[END_REF][START_REF] Zhang | Total-evidence dating under the fossilized birth-death process[END_REF]. This leads to a statistically coherent prior on divergence times, where the variance associated with node ages reflects the incompleteness of the fossil record, as well as uncertainty associated with the placement of fossil taxa in the phylogeny. Bayesian inference using the FBD process as a tree prior also allows for reliable estimation of the diversification and sampling parameters, although at least one of the parameters is required to be fixed by the user for the model to be identifiable [START_REF] Stadler | On incomplete sampling under birth-death models and connections to the samplingbased coalescent[END_REF][START_REF] Gavryushkina | Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration[END_REF]. In most empirical applications, the extant sampling proportion is the best known and is thus chosen to be fixed. The FBD model has been successfully applied to datasets of living and fossil taxa [START_REF] Schuster | Divergence times in demosponges (Porifera): first insights from new mitogenomes and the inclusion of fossils in a birth-death clock model[END_REF][START_REF] Šmíd | Calibrating the tree of vipers under the fossilized birth-death model[END_REF][START_REF] Daniel | Ancient crested penguin constrains timing of recruitment into seabird hotspot[END_REF][START_REF] Pohle | Early cephalopod evolution clarified through bayesian phylogenetic inference[END_REF].

The initial implementation of the FBD model assumed constant diversification (birth and death) and fossil recovery rates through the entire phylogeny [START_REF] Heath | The fossilized birth-death process for coherent calibration of divergence-time estimates[END_REF]. However, a wide range of factors, from biological and geological processes to collection practices, contribute to the probability a given organism will be sampled in the fossil record [START_REF] Kidwell | The quality of the fossil record: Implications for evolutionary analyses[END_REF][START_REF] Smith | The ties linking rock and fossil records and why they are important for palaeobiodiversity studies[END_REF][START_REF] Roger Bj Benson | Biodiversity across space and time in the fossil record[END_REF][START_REF] Nussaïbah B Raja | Colonial history and global economics distort our understanding of deep-time biodiversity[END_REF]. As a consequence, fossil recovery potential varies substantially across time, space and taxa. Other variables, such as environmental conditions and phenotypic traits, contribute to variation in diversification (speciation and extinction) rates.

Extensions of the FBD process have integrated these variations into the model [START_REF] Gavryushkina | Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration[END_REF][START_REF] Zhang | Total-evidence dating under the fossilized birth-death process[END_REF], however, many empirical analyses are still done under the constantrate assumption, for several reasons such as computational cost, lack of precise knowledge of rate changes, and ease of setup and interpretation.

One challenge in integrating fossil specimens in FBD analyses is that genetic information degrades with time. Thus, molecular sequences are generally unavailable for fossils, unlike extant species. As a result, additional information needs to be added to the inference to allow fossils to be placed in the tree topology. In previous research, this information has been provided in two different forms, which can be used separately or in combination. The first is topological constraints, which are constraints added by the user specifying that certain subsets of tips, extant or extinct, need to be monophyletic clades in the inferred phylogeny. These constraints are generally based on the taxonomy, with the constrained clades corresponding to genera and/or higher classifications. Topological constraints use information which is usually already available, and do not add computational complexity or cost to the inference. However, they do not easily accommodate uncertainty in the taxonomy, which is present even for well-known crown groups [START_REF] Marx | Cetacean Paleobiology[END_REF]. Another approach is so-called total-evidence analyses [START_REF] Ronquist | A total-evidence approach to dating with fossils, applied to the early radiation of the Hymenoptera[END_REF], which uses a morphological data matrix, with data for both fossil and extant specimens, in addition to the molecular alignment. While morphological characters can be discrete or continuous, matrices used in total-evidence approaches are usually composed exclusively of discrete (or discretized) characters. The morphological matrix is added to the inference along with a morphological substitution model and a morphological clock model, and contributes to the phylogenetic likelihood. Although total-evidence approaches better account for the underlying empirical data, they are more costly both in the time and effort required to assemble the matrix and in the added computational cost of the inference. In addition, the accuracy and precision of the inference has been found to be strongly dependent on the size of the morphological matrix [START_REF] Barido-Sottani | Ignoring fossil age uncertainty leads to inaccurate topology and divergence time estimates in time calibrated tree inference[END_REF].

Although these two approaches can in theory be combined in the same analysis, in practice they are often viewed as separate alternatives (for instance, [START_REF] Šmíd | Calibrating the tree of vipers under the fossilized birth-death model[END_REF] uses constraints, while Thomas et al. ( 2020) uses a total-evidence approach). However, which method to choose when both are available remains unknown. In addition, we expect that the presence of uncertainty in the information used to place fossils (in the form of imprecise constraints or small morphological matrices) will lead to lower accuracy, but this has not been precisely evaluated. In this work, we compare and evaluate different approaches to add fossil placement information in FBD analyses, both in ideal conditions and in datasets including uncertainty or even errors. We also evaluate the performance of combining both topological constraints and morphological data in the same inference. Finally, since the phylogeny estimation in a Bayesian inference is impacted both by the data and by the chosen tree model, we also explore the impact of model misspecification on the accuracy of inferences. By using datasets containing variation in fossil recovery rates or speciation rates, we study whether some fossil placement strategies are more robust to unaccounted-for rate variations in the underlying processes.

Methods

Simulations

Trees and fossils

Our goal was to assess the performance of the FBD process, examining the impact of phylogenetic uncertainty and model misspecification. The parameters of the simulation were constrained to reflect values obtained for marine invertebrates.

Trees and fossils were simulated using the R packages TreeSim [START_REF] Stadler | Simulating trees with a fixed number of extant species[END_REF]FossilSim (Barido-Sottani et al., 2019b). The simulations were conditioned on the number of extant tips (n = 25), with constant speciation rate (λ = 0.11), extinction rate (µ = 0.1) and fossil recovery rate (ψ = 0.03). We assumed complete sampling at the present (i.e., the probability of extant species sampling ρ = 1). Parameters were selected such that the expected origin time was 250 Myr and the expected number of fossils was 100. Simulated data sets were filtered to select for trees with an origin time between 225 and 275 Myr and with between 80 and 120 fossils, and simulations which did not fit these two criteria were discarded. In addition, we simulated two datasets with incomplete extant sampling, presented in more details in the Supplementary Materials.

We also simulated two sets of trees integrating two different sources of rate variation. Set (a) used constant speciation and extinction rates (λ = 0.11, µ = 0.1) and variable fossil recovery rates (ψ 1 = 0.02 and ψ 2 = 0.04), while set (b) used variable speciation rates (λ 1 = 0.12 and λ 2 = 0.06), a constant extinction rate (µ = 0.1) and variable fossil recovery rates (ψ 1 = 0.02 and ψ 2 = 0.04).

In both sets, the non-variable parameters (i.e., the speciation rate in set (a) and the extinction rate in sets (a) and (b)) were fixed to the same value as the constant birth-death simulations. In both variable sets, the variation in rates was linked to a trait ch that was simulated along each tree under a Brownian motion process. Traits were simulated with a variance σ = 0.01 in both sets, and with an initial value of ch origin = 2.25 for set (a) and ch origin = 1.5 for set (b). Trait values were then assigned to two discrete types: values ch < 0 were assigned to type 1 (with rates ψ 1 and λ 1 ) and values ch > 0 were assigned to type 2 (with rates ψ 2 and λ 2 ). The final trees were kept if each trait value was assigned to at least 10% of the tips. The parameters of the BM process were calibrated such that the final sampled trees contained an average of 8 to 11 trait changes.

We generated 50 replicate phylogenies for each of the three simulation conditions, resulting in 150 trees in total.

Characters

We simulated both molecular and morphological alignments for all 50 phylogenies in each set of trees. Molecular sequence alignments of 1000 sites were simulated under the HKY +Γ model with five discrete gamma rate categories (α = 0.25), using the implementation of Seq-Gen available in the R package phyclust. Branch rates were simulated under a lognormal uncorrelated clock model. For each tree replicate the average substitution rate was sampled from a gamma distribution with an expected value of 0.01 substitution per My, and shape and scale parameters = 2 and 0.5, respectively. The log of this rate was then used to define the mean of a lognormal distribution with variance = 0.01 from which branch specific rates were independently drawn.

Based on the results shown in [START_REF] Barido-Sottani | Ignoring fossil age uncertainty leads to inaccurate topology and divergence time estimates in time calibrated tree inference[END_REF], which showed that total-evidence approaches using only 30 characters led to unreliable estimates, we selected alignment sizes of 50 and 300 characters for our morphological data matrices, corresponding to low-sampling and highsampling scenarios, respectively. These matrices were simulated under a binary state Mk model with a strict clock and a rate of 0.1 substitutions per My, using the sim.char function provided by the R package geiger.

Inference

We analysed each replicate using Bayesian phylogenetic inference in the BEAST2 framework [START_REF] Bouckaert | BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis[END_REF], under the constant rate FBD process implemented in the package Sampled Ancestors (SA) [START_REF] Gavryushkina | Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration[END_REF].

We examined the impact of eight different ways of incorporating the phylogenetic uncertainty associated with fossils:

1. Fossil samples without character data were assigned to the nearest extant ancestral node using clade constraints (designated as "correct constraints")(Fig. 1C). (1) but with 2% of fossils assigned to the wrong node, selected at random (designated as "constraints with errors")(Fig. 1D).

As in

3. Fossil samples without character data were assigned to the node above the nearest ancestral node (designated as "imprecise constraints")(Fig. 1E).

4. As in (1) but 5 nodes picked at random in the tree were collapsed, meaning all constraints below that node were removed (designated as "collapsed constraints")(Fig. 1F).

5. Character data (= 50 characters) was included for both fossil and extant samples, with no additional constraints (designated as "total-evidence with n =50")(Fig. 1G).

6. Character data (= 300 characters) was included for both fossil and extant samples, with no additional constraints (designated as "total-evidence with n =300")(Fig. 1G).

7. Both character data (= 300 characters, as in ( 6)) and constraints with errors (as in (2)) were applied (designated as "combined evidence with errors on constraints") 8. Both character data (= 300 characters, as in ( 6)) and imprecise constraints (as in (3)) were applied (designated as "combined evidence with imprecise constraints") Configuration ( 1) corresponds to a situation with perfect information on fossil taxonomy (Fig. 1C), while configuration (2) includes a low amount of misplaced fossils (Fig. 1D). Configurations ( 3) and ( 4) include uncertainty in fossil placement in two different ways: in configuration (3) all clades have some amount of uncertainty in fossil assignments (Fig. 1E), in configuration (4) a few clades have no information on fossil placement within the clade, while the other clades are known perfectly (Fig. 1F). Configurations ( 5) and ( 6) are total-evidence analyses, with respectively a low or high number of morphological characters for each fossil (Fig. 1G). Finally, configurations ( 7) and ( 8) represent a situation where morphological data is used to complement incomplete or uncertain taxonomy information. These last two configurations were only used in the dataset with no rate variation.

The FBD process was parameterized using the 'canonical' (speciation, extinction, sampling) parameterization, with an exponential prior (mean = 1.0) on the birth, death, and sampling parameters.

A uniform prior U (0, 1000) was used for the origin time of the process. We deliberately chose uninformative priors for these parameters in order to minimize the influence of the priors on the results. The extant sampling proportion was fixed to the true value of ρ = 1. A lognormal prior was used for the HKY parameter κ (mean = 1 and sd = 1.25), with a uniform prior on the state frequencies. An exponential prior was used for the alpha parameter of gamma distributed rates (mean = 1). The Mk model was applied to the binary character data. An exponential prior (mean = 1) was used for the mean of the uncorrelated relaxed clock model and a gamma prior (α = 0.5396 and β = 0.3819) on the standard deviation.

Analyses were run for a minimum of 200,000,000 generations, sampling every 10,000 generations and discarding 10% as burnin. We assessed convergence by calculating ESS values for the likelihood, prior and posterior using the R package coda [START_REF] Plummer | Coda: convergence diagnosis and output analysis for mcmc[END_REF]. Since we were interested in the accuracy of the topology, we also calculated the median ESS values for phylogenetic splits in the posterior. This was done using the R package convenience [START_REF] Guimaraes | Convergence assessment for bayesian phylogenetic analysis using mcmc simulation[END_REF] based on a subsample of 2000 posterior trees (it was necessary to downsample for computational reasons). Finally, we calculated an approximate ESS for the overall topology using the R package rwty [START_REF] Dan L Warren | Rwty (r we there yet): an r package for examining convergence of bayesian phylogenetic analyses[END_REF]. We aimed to achieve ESS values > 200. This ensures that the standard error of the estimate is < 2% (Fabreti and Höhna, 2022, see table 1), which we regard as efficient for simulations. If ESS values were << 200, analyses were longer (this was only necessary for 23 replicates). 13 analyses that failed to converge after 2,000,000,000 generations were excluded.

We note that nine of these are analyses in which fossils have no character data and some nodes are collapsed, with fossils assigned to the next nearest ancestor (Configuration 4, Fig. 1F). The maximum number excluded for a given simulation and inference scenario was 3 (out of 50 replicates for each configuration). ESS values are provided in a table in the supplementary material.

For each replicate, we calculated the relative error of divergence time estimates (as the absolute difference between the median estimate and the true value, divided by the true value), averaged over all extant nodes, and the 95% HPD coverage, i.e., the proportion of extant nodes for which the true age was contained within the 95% Highest Posterior Density interval. We also measured the agreement between the estimated phylogeny and the truth using two separate metrics, the normalised Robinson-Foulds distance, calculated using the R package phangorn and the quartet similarity, calculated using the R package Quartet. Both metrics were averaged over all posterior samples, for both the full tree including fossils and the reconstructed extant phylogeny.

We also performed a separate unconstrained (i.e., non-FBD) tree inference using RevBayes. Since these analyses contain no age information the trees produced are non-time calibrated and have branch lengths in substitutions per site rather than in units of time. We used a uniform tree prior on the topology, with an exponential prior (λ = 10, mean = 0.1) on the branch lengths. We used the same settings as above for substitution models. We ran 4 independent chains for 200,000 generations, sampling every 200 generations, discarding 10% as burnin and combining the output.

We assessed convergence by calculating the ESS values as above and examined the variance across chains. All runs converged. We measured the normalised Robinson-Foulds distance between the estimated phylogeny and the truth, averaged over all posterior samples, and compared those to the results obtained using the FBD inference.

Results

The results on the datasets simulated under constant rates, variable fossilization rates, and variable birth and fossilization rates are shown in Figures 2 and3 Under a constant-rate process, we can see that the extant topology is well recovered under all methods of fossil placement. Divergence times are similarly well recovered across all methods, with the marked exception of constraints which contain errors. Indeed, the relative error is much higher and the coverage much lower in the presence of errors. Combining morphological data and constraints improves the divergence time estimates markedly compared to the same constraints alone. In particular, our results show that adding a morphological alignment can compensate somewhat for the presence of errors in the topological constraints, although the accuracy remains much lower than in any of the configurations without errors.

As expected, fossil placement methods also impact the accuracy of the estimate of the full phylogeny. The best performance when using only one method of fossil placement is obtained with correct constraints, followed by imprecise constraints and total-evidence with a high number of morphological characters, whereas the worst estimates are obtained using total-evidence with a low number of morphological characters. Consistent with the results on the divergence time estimates, combining morphological data with constraints significantly improves the accuracy of the phylogeny estimate. Combined inferences perform better at reconstructing the full phylogeny than either of the fossil placement strategies alone, including when using complete correct constraints.

We see very similar patterns in datasets which include rate variation. One notable difference is that relative errors in extant divergence times are higher the more variation is included in the dataset for all approaches using topological constraints, and particularly for constraints with errors.

Interestingly, total-evidence approaches show the same levels of relative error in extant divergence times regardless of rate variation in the simulation models. One possible explanation is that totalevidence approaches include more uncertainty than topological constraints and so are able to better compensate for model mismatches between the simulation and the inference. The effect of simulated rate variation is almost null on the other metrics, i.e., the coverage of extant divergence times and the RF distances for extant and full phylogenies.

Figure . 4 shows the results obtained using the FBD versus non-FBD (i.e., without fossil age information) inference. The comparison shows that trees recovered under the FBD model are more accurate, in terms of Robinson-Foulds distance -this applies to both the extant topology estimated using the molecular data only or the total-evidence matrix and the full topology. This result holds even with the use of erroneous fossil constraints and model misspecification under the FBD. Interestingly, the unconstrained topology estimates are more accurate for trees simulated using constant diversification and sampling rates than trees simulated with rate variation. This might reflect a better correspondence between the uniform tree model and a constant birth-death sampling pro-cess. The inclusion of fossil taxa increases the accuracy of the extant tree. More morphological 260 character data also leads to improved accuracy of both the extant and full trees. We then set up clade constraints according to different rules: in (C), all clade constraints are correct and complete. In (D), some fossils are assigned to the wrong clade. In (E), all clades have a low level of taxonomic uncertainty. In (F), some subclades are fully unknown while the rest of the tree is fully known. Finally, we also perform two total-evidence analyses (G). , with no fossil age information) for simulations under a constant birth-death process (blue), variation in fossilization rates (pink) or variation in birth and fossilization rates (green). Average Robinson-Foulds distance between posterior samples and the true tree, for the extant tree only (A) or the full tree including fossil samples (B). All measures show the average and standard deviation across all replicates for different fossil placement approaches.

Discussion

The fossilised birth-death process offers flexible opportunities to reconstruct dated phylogenies incorporating fossils under a mechanistic framework. Here, we explored different options for handling the phylogenetic or taxonomic uncertainty associated with fossil samples using simulations.

We found that alternative approaches to including fossils had relatively little impact on the accuracy of the extant topology (based on RF distances, Fig. 2C-D). Overall, these results are quite consistent with expectations: in FBD inferences, divergence times are mostly informed by fossil ages, so variations in the position of fossils will strongly impact these estimates. On the other hand, the placement of extant taxa in the phylogeny is driven primarily by the molecular sequence alignment, and thus is mostly independent from the handling of fossil specimens. Using fixed constraints to correctly assign fossils to nodes produced more accurate trees, compared to trees recovered using morphological characters to place the fossils. However, this represents an idealised scenario, as reliable and precise taxonomic information used to inform constraints will often not be available (see below for further discussion on data quality and availability). The results obtained using nontime calibrated tree inference show that trees recovered with the inclusion of fossils with character data are more accurate than trees based on the extant taxa only (Fig. 4A-B). We also confirm that results are even better when the fossil age information is taken into account using the FBD process.

Previous simulations studies have also shown that the inclusion of fossils and age information can improve the accuracy of the topology among extant taxa (Mongiardino [START_REF] Nicolás | Fossils improve phylogenetic analyses of morphological characters[END_REF].

Different approaches to handling taxonomic uncertainty had a more notable impact on divergence times. We found that the best possible scenario, in terms of recovering accurate node ages (relative error and coverage), was to assign fossils to the correct nearest node in the extant tree (Fig. 2A-B). This is unsurprising, as this is equivalent to fixing large portions of the topology using correct monophyletic clade constraints. But we further show that the use of less precise clade constraints (i.e., assigning fossils to larger or more inclusive clades, Fig. 1E) also recovers accurate estimates, similar to those obtained inferring the position of fossils based on morphological character data.

Similarly, [START_REF] Heath | The fossilized birth-death process for coherent calibration of divergence-time estimates[END_REF] found using simulations that less precise constraints increased the variance but did not reduce accuracy of posterior divergence estimates (see also O'Reilly and Donoghue ( 2020)). However, our results show that when only a small proportion (2%) of fossils are assigned to incorrect nodes, overall accuracy decreases. Assigning larger proportions of fossils to incorrect nodes will reduce accuracy further. Preliminary runs using 10% of fossils assigned to incorrect nodes lead to coverage < 0.1 (results not shown). Together, our results suggest that uncertainty is a much less critical issue for inference than outright errors. Thus, we recommend that when in doubt, it is better to err on the side of caution and to either use larger and more inclusive clade constraints, since this does not compromise the overall accuracy of results and only leads to a small overall decrease in precision. This will be especially important to consider for taxa whose position within the stem vs crown group is disputed, since early crown group representatives can be difficult to distinguish from stem group taxa [START_REF] Robert | Bias and sensitivity in the placement of fossil taxa resulting from interpretations of missing data[END_REF]. If in any doubt, we suggest using 'total-group' constraints, rather than assigning fossils to the stem or crown group a priori.

Alternatively, we can infer the position of fossils using morphological data when available. This is preferential since the posterior will best reflect uncertainty associated with the placement of fossil taxa. In addition, if clade constraints are overly conservative (i.e., based on very large clades or higher level taxonomic divisions) there will be too much uncertainty and the analyses might fail to converge. 'Imprecise constraints' was the simulation scenario in which we had to exclude the largest number of replicates (= 4/50) due to lack of convergence.

We also examined the impact of model violation on the accuracy of topology and divergence times recovered under the FBD model, by varying sampling rate only or both sampling and diversification rates. Not accounting for rate variation in these parameters has a modest impact on the accuracy of the topology -again reflecting the fact that the signal for topology predominantly comes from the character data (Fig. 2C-D). Model misspecification has a much more discernible impact on the divergence estimates (Fig. 2A-B). This is because the signal for divergence times largely comes from the distribution of fossil sampling times and relies more on the birth-death sampling model being correct. However, in most scenarios coverage for node ages remains high. These general results also match the findings of previous simulation studies that explored model misspecification under the FBD process [START_REF] Heath | The fossilized birth-death process for coherent calibration of divergence-time estimates[END_REF]. There are challenging identifiability issues associated with birth-death processes [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF][START_REF] Louca | Fundamental identifiability limits in molecular epidemiology[END_REF]. In particular, [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF] identified 'congruence classes', within which infinitely many diversification histories can have the same likelihood. This means that applying oversimplified birth-death models can result in highly misleading results. Our findings show that reliable phylogenies can be obtained under the FBD process, even if the underlying process is more complex than the model used for inference (see also [START_REF] Koch | Inaccurate fossil placement does not compromise tip-dated divergence times[END_REF]). The most notable exception are the results obtained with erroneous fossil constraints -when a portion of fossils are incorrectly placed in the tree, the relative negative impact of model misspecification is worse. Since we can rarely be certain of the underlying diversification or sampling process, this further underscores the need to use taxonomic constraints with extreme caution. Furthermore, our results demonstrate that total-evidence analyses are more robust against model misspecification, adding to the list of benefits associated with this approach.

Ordinarily clade constraints and 'total-evidence' analyses are treated as two distinct approaches, but they can easily be combined. Several previous simulation studies have shown that the accuracy of parameters estimated under the FBD process increases with fossil sampling [START_REF] Heath | The fossilized birth-death process for coherent calibration of divergence-time estimates[END_REF][START_REF] Barido-Sottani | Ignoring fossil age uncertainty leads to inaccurate topology and divergence time estimates in time calibrated tree inference[END_REF], and our results confirm that the combined approach yields more accurate estimates than either fossil placement strategy used alone. Combined inferences perform particularly well when reconstructing the phylogeny itself, and the addition of morphological data can compensate partly for the presence of errors in the constraints when estimating the divergence times. The combined strategy seems particularly relevant for empirical analyses. Indeed, in many datasets, we might only have abundant morphological character data associated with a subset of fossils. For example, we might have fossils with fewer traits (e.g., fossil cephalopod shells), alongside rarer but exceptionally preserved specimens with more informative traits (e.g., fossil cephalopods associated with jaws and or soft-tissues). Shell characteristics that are diagnostic of higher-level taxonomic groups could be used to define clade constraints for specimens associated with a smaller amount of data. Similarly, in palaeobotany, rarely preserved tissues (e.g., flowers or seeds) can have a small number of traits (i.e., 2-3) that would be considered too few to construct a matrix for phylogenetic analysis, but that are nonetheless considered definitive synapomorphies of certain plant groups. Another example could be ammonoid aptychi or lower jaws, which have been mapped to phylogenies based on a small number of specimens with known jaws, relative to those where shells/moulds are known [START_REF] Engeser | Phylogeny of the aptychi-possessing neoammonoidea (aptychophora nov., cephalopoda)[END_REF]. Parasites might be an especially good example in this context. Many groups of parasites are highly diverse today but most have poor preservation in the past, associated with two distinct fossil datasets -rare exceptionally preserved specimens that can be used to build a morphological matrix versus more common eggs or trace fossils that can only be assigned to higher groupings [START_REF] De Baets | The fossil record of parasitism: Its extent and taphonomic constraints[END_REF], potentially better suited to the use of clade constraints [START_REF] Warnock | The molecular clock as a tool for understanding hostparasite evolution[END_REF].

We note that our simulation design represents an idealised scenario in several ways. For instance, we only simulated binary characters for our morphological alignments, whereas empirical data matrices often include multi-state characters or characters with hierarchical state dependencies. In addition, our matrices had no missing data, which is unrealistic for empirical datasets including fossils. These two factors likely contribute to the good overall performance of the total-evidence analyses in our study. However, similar results were found recently by [START_REF] Koch | Inaccurate fossil placement does not compromise tip-dated divergence times[END_REF] in a totally different set of simulations that employed a more complex model of morphological evolution (linked to the diversification process) and incorporated missing data. Both studies show that inaccurate or imprecise estimated fossil placement does not have an adverse impact on divergence times, as opposed to the use of inaccurate constraints. Finally, our main results were obtained assuming complete extant sampling, which is unrealistic in many empirical datasets, especially for larger clades [START_REF] Matschiner | Selective sampling of species and fossils influences age estimates under the fossilized birth-death model[END_REF]. Our results on incompletely sampled datasets, presented in the Supplementary Materials, show that divergence time estimates are largely unaffected by incomplete sampling. Phylogeny estimates show a noticeable decrease in accuracy when incompletely sampled datasets are combined with a total-evidence approach, but not when run with complete topological constraints. A more comprehensive exploration of the impact of incomplete sampling, particularly in combination with rate variations in the diversification and fossilization processes, remains to be done.

Here, we considered a scenario in which we have an extant clade. For fossil only datasets, typically we use the FBD model with morphological data only and no topological constraints. Previous simulations have shown that trees generated with matrices that are typical of many fossil groups (i.e., 30 characters) will be highly uncertain [START_REF] Barido-Sottani | Ignoring fossil age uncertainty leads to inaccurate topology and divergence time estimates in time calibrated tree inference[END_REF]. Clade constraints could also be useful in this context, provided these can be defined with a high degree of confidence. As an example, a backbone constraint based on the analysis of molecular data could be applied to a tree that consists mostly of fossils or is based on morphology, but where the morphology never recovers the established molecular phylogeny. For example, cyclostomes [START_REF] Near | Conflict and resolution between phylogenies inferred from molecular and phenotypic data sets for hagfish, lampreys, and gnathostomes[END_REF] or many mammalian superorders [START_REF] Foley | Mammal madness: is the mammal tree of life not yet resolved?[END_REF][START_REF] Robin | Improvements in the fossil record may largely resolve current conflicts between morphological and molecular estimates of mammal phylogeny[END_REF] are never recovered as monophyletic on the basis of morphology alone. In other cases, we might want to look at a large clade that contains too many fossil species for all of them to be included in a single analyses, but where previous studies show unequivocal support for certain subclades. In this scenario, clade constraints could help stabilising the topology, without the need to collect additional morphological data while remaining computationally efficient.

Further studies are needed to establish the best fossil placement strategies in situations where the classical FBD approach might not be suitable. One such situation is when different types of evidence are associated with the presence of a group. For instance, worm eggs are sometimes assigned to higher taxonomic ranks based on the presence of specific structures [START_REF] Hugot | Discovery of a 240 million year old nematode parasite egg in a cynodont coprolite sheds light on the early origin of pinworms in vertebrates[END_REF], but these traits can have little to do with the adult morphology of exceptionally preserved worms. Other examples include trace fossils (e.g., where they predate body fossils), molecular fossils (e.g., biomarkers), or exceptional preservation of non-adult ontogenetic stages (e.g., small larvae in early phosphatic or Orsten-type preservation, see [START_REF] Maas | The 'orsten'-more than a cambrian konservat-lagerstätte yielding exceptional preservation[END_REF]). These alternative lines of evidence cannot be reconciled easily within an analysis using a single morphological matrix, even in combination with clade constraints. One solution in this scenario could be to partition different types of fossil characters and use a multi-loci approach such as the multi-species coalescent with fossils (FBD-MSC) model [START_REF] Huw | Novel Integrative Modeling of Molecules and Morphology across Evolutionary Timescales[END_REF].

For some groups (e.g., soft-bodied worms or unicellular organisms), it might not be possible or practical to collect more data or obtain reliable information about broader level taxonomy, due to both incomplete preservation and/or the amount of work associated with collecting morphological data. In this situation, extended versions of the FBD process can help make the best use of the available data. For instance, the occurrence birth-death process allows fossils associated with morphological character data and fossils associated with age data only to be assigned different sampling rates [START_REF] Andréoletti | The occurrence birth-death process for combined-evidence analysis in macroevolution and epidemiology[END_REF]. Under this model, occurrence data do not need to be constrained to any part of the tree, but can improve overall FBD parameter estimates.

Beyond inferring dated trees, Soul and Friedman (2015) showed that trees constructed using higher level taxonomy and dated using time-scaling methods can be successfully used for phylogenetic comparative analyses. Although we emphasise the need to take a very cautious approach when using clade constraints, the FBD model could also be used in this context, with the added benefit that the output better reflects uncertainty in fossil ages, node ages and phylogenetic uncertainty.

Irrespective of the approach to handling taxonomic or phylogenetic uncertainty using the FBD process, core taxonomic research, including morphological data collection and species descriptions, remains essential. The issues we observed with the erroneous placement of fossils mirror those identified previously in node dating studies, where node calibrations based on inaccurate fossil as- et al., 2012). Authors are recommended to provide an explicit set of statements justifying the assignment of a fossil to a given node, with reference to up-to-date phylogenetic analyses incorporating the relevant taxa or to a set of unequivocal synapomorphies. This rigorous and transparent approach to defining clade constraints is directly applicable to analyses that employ the FBD process, along with the criteria used to justify fossil ages. How age uncertainty is handled in analyses using the FBD process also has important implications for the accuracy of both divergence times and topology (Barido-Sottani et al., 2019a, 2020). Developments of phylogenetic models used in palaeobiology are no substitute for the expertise contributed by fundamental systematics, taxonomic and stratigraphic research. Our results reiterate the need for increased and direct support for taxonomy-based projects [START_REF] Britz | Funding, training, permits-the three big challenges of taxonomy[END_REF][START_REF] Michael | The taxonomic impediment: a shortage of taxonomists, not the lack of technical approaches[END_REF]. Improving approaches to phylogenetic dating requires both advanced methodological and empirical perspectives.

  for the divergence time estimates and the Robinson-Foulds distances, and Figures S1 and S2 for the quartet similarities. The Robinson-Foulds distance and the quartet similarity led to very similar patterns.

Figure 1 :

 1 Figure 1: Schematic representation of the different analyses performed in this study. First, we simulate a full tree with fossil samples (A), from which the true extant tree can be obtained (B).

Figure 2 :Figure 3 :Figure 4 :

 234 Figure2: Results of FBD analyses for simulations under a constant birth-death process (blue), variation in fossilization rates (pink) or variation in birth and fossilization rates (green). Absolute relative error (A) and 95% HPD coverage (B) of divergence time estimates, averaged over all nodes in the extant phylogeny. Average Robinson-Foulds distance between posterior samples and the true tree, for the extant tree only (C) or the full tree including fossil samples (D). All measures show the average and standard deviation across all replicates for different fossil placement approaches.

  signments have been shown to result in large errors in divergence time estimation(e.g., Phillips et al., 2009). These issues spurred the development of 'best practices' for justifying fossil calibrations for node dating, taking into account both phylogenetic and fossil age uncertainty (Parham

6 Acknowledgments JBS was supported by funds from the National Science Foundation (USA), grants DBI-1759909 and DEB-1556615 and from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Sk lodowska-Curie grant agreement No. 101022928. KDB was supported by the I.3.4 Action of the Excellence Initiative -Research University Programme at the University of Warsaw (Project: PARADIVE).

Author contributions

JBS and RCMW conceived of the study and performed the analyses. All authors contributed to the design of the study, the interpretation of results and the writing of the manuscript.

Data availability

The code used to perform the simulations and analyses, as well as the simulated datasets are available on Zenodo (DOI: 10.5281/zenodo.6884064).