# Exploring congruent diversification histories with flexibility and parsimony

Jérémy Andréoletti, Hélène Morlon

► To cite this version:

Jérémy Andréoletti, Hélène Morlon. Exploring congruent diversification histories with flexibility and parsimony. Methods in Ecology and Evolution, 2023, 14 (12), pp.2931-2941. 10.1101/2023.07.26.550618. http://hal-04173637

HAL Id: hal-04173637

https://hal.science/hal-04173637v1

Submitted on 13 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# EXPLORING CONGRUENT DIVERSIFICATION HISTORIES WITH
FLEXIBILITY AND PARSIMONY

Jérémy Andréoletti

Institut de Biologie - École Normale Supérieure - Université PSL

j Jeremy.andreoletti@gmail.com

Hélène Morlon

Institut de Biologie - École Normale Supérieure

Université PSL, CNRS, INSERM

## ABSTRACT

1. Using phylogenies of present-day species to estimate diversification rate trajectories – speciation and extinction rates over time – is a challenging task due to non-identifiability issues. Given a phylogeny, there exists an infinite set of trajectories that result in the same likelihood; this set has been coined a congruence class. Previous work has developed approaches for sampling trajectories within a given congruence class, and suggested that rapid changes in speciation or extinction rates are conserved across the class.

2. We introduce a new method for exploring congruence classes, that we implement in the R package CRABS. Whereas existing methods constrain either the speciation rate or the extinction rate trajectory, ours provides more flexibility by sampling speciation and extinction rate trajectories simultaneously. We also implement a filtering step that allows selecting the most parsimonious trajectories within a class.

3. We demonstrate the utility of our new sampling strategy using a simulated scenario. Next, we apply our approach to the study of mammalian diversification history. We show that rapid changes in speciation and extinction rates need not be conserved across a congruence class, but that selecting the most parsimonious trajectories shrinks the class to concordant scenarios.

4. Our approach opens new avenues both to truly explore the myriad of potential diversification histories consistent with a given phylogeny, embracing the uncertainty inherent to phylogenetic diversification models, and to select among these different histories. This should help refining our inference of diversification trajectories from extant data.
## 1 Introduction

The study of diversification rates, or the rates at which species arise and go extinct over time, has long been a central topic in evolutionary biology. Early approaches for inferring diversification rates relied on simple heuristic methods, based primarily on observed species ranges in the fossil record (Raup, 1978; Foote, 1994; Jablonski & Chaloner, 1994), and have gradually been expanded to fit into more robust statistical frameworks (Foote, 2000; Alroy, 2014; Silvestro et al., 2014). Another major innovation was to incorporate phylogenetic information from molecular data (Nee et al., 1994) and to model diversification as a continuous branching process of speciation, extinction and sampling (Stadler, 2010). These mathematical developments paved the way for a series of likelihood-based approaches to estimating diversification rates with increasingly complex models, taking into account both rate variation over time (Morlon et al., 2011; Stadler, 2011) and rate heterogeneity across lineages (Alfaro et al., 2009; Rabosky, 2014; Maliet et al., 2019; Barido-Sottani et al., 2020).

However, in recent years, a fundamental problem in statistical inference has been brought to the forefront of diversification rate estimation: the lack of parameter identifiability. For some birth-death diversification models, it is not possible to uniquely determine the underlying parameters of the diversification process from observed data. This issue had already been raised by Kubo & Iwasa (1995), Stadler (2009) and Lambert & Stadler (2013), but Louca & Pennell (2020) demonstrated in full generality that, for a homogeneous birth-death process with speciation and extinction rates that vary freely in time, an infinite combination of speciation and extinction rate curves (usually noted $\lambda(t)$ and $\mu(t)$ respectively) yield the exact same likelihood for a given phylogeny of extant species, even in the limit of infinite data. They refer to the infinite set of trajectories $(\lambda(t), \mu(t))$ that are strictly indistinguishable from the phylogeny alone as a "congruence class". Louca et al. (2021) later extended this non-identifiability result to birth-death processes that also include a fossil sampling rate.

Three main avenues have been followed to address this finding. First, Louca & Pennell (2020) proposed to restrict inferences to identifiable reparametrizations of the model (see also Stadler (2013)), though these "pulled rates" can be challenging to interpret (Helmstetter et al., 2022). The second approach is to work with hypotheses-driven models that are not as flexible as those used by Louca & Pennell (2020), or to add regularity constraints (e.g. priors in a Bayesian inference context) that make the process identifiable (Morlon et al., 2022). Many models that are not as flexible as those used by Louca & Pennell (2020) are indeed identifiable (Legried & Terhorst, 2022a,b). The final option is to embrace the non-identifiability results and to develop methods to explore scenarios within a given congruence class, in order to assess which temporal patterns are conserved across the congruence class and which are not (Höhna et al., 2022; Kopperud et al., 2023). The latter avenue is the one we investigate here, and we link it to the second approach.

Current approaches for exploring congruent classes are based on results from Louca & Pennell (2020), which showed that each congruence class can be uniquely defined by the speciation rate at present $\lambda_0$ and a curve called the "pulled diversification rate" curve, defined as $r_p(t) := \lambda(t) - \mu(t) + \frac{1}{\lambda(t)} \frac{d\lambda(t)}{dt}$. Given $\lambda_0$ and $r_p$, one can explore trajectories in the congruence class by
choosing any trajectory for either the speciation rate or the extinction rate and finding the trajectory for the other rate using this equation. In practice, this can be done with a function called `congruent_hbds_model` from the R package `castor` (Louca & Doebeli, 2018). Höhna et al. (2022) developed this approach further with two methods for sampling nearly congruent trajectories based on the following linear approximation (for the time step $t_i$ at iteration $i$): $r_p(t_i) = \lambda_i - \mu_i + \frac{1}{\lambda_i} \frac{\lambda_i - \lambda_{i-1}}{\Delta t}$. For simplicity we will refer to "congruent trajectories" in what follows, although it should be kept in mind that the trajectories are not exactly congruent, given the linear approximation. The first method is similar to `congruent_hbds_model`, as it involves choosing a trajectory for either the speciation or the extinction rate, and then computing the other trajectory. The second method consists in sampling $\lambda(t)$ or $\mu(t)$ more freely in a series of specified function classes, e.g. constant extinction rates with various initial $\mu_0$ values ($\lambda_0$ is invariant in each congruence class), Gaussian or Horseshoe Markov random fields (see Magee et al. (2020)), or even a purely uniform sampling at each time step. This aims at truly sampling scenarios in the congruent class, as opposed to hand-picking few of them. The authors implemented these approaches in a new R package called CRABS (Congruent Rate Analyses in Birth-Death Scenarios). In a follow-up paper, Kopperud et al. (2023) used CRABS to investigate the speciation curves resulting from sampling extinction curves, as sampling speciation curves to recover congruent extinction curves is trickier (see below). With this approach, the authors showed that moderate rate changes are not necessarily recovered across congruent models, but that sharp shifts often are.

One main limitation of these approaches is that they strongly constrain the rate trajectory that is fixed or sampled ($\lambda$ or $\mu$), as opposed to the trajectory that is recovered through the congruence class, hence limiting the diversity of patterns that can be observed. Even the CRABS option of choosing a series of uniformly sampled points imposes the absence of any recurrent pattern on the rate trajectory that is sampled. As we illustrate in the paper, this can be problematic. First, if only one strategy is used (for example sampling mostly extinction rates as in Kopperud et al. (2023)), we may misleadingly conclude that the inference is robust across the congruence class. If instead several alternative sampling strategies are used, as recommended by Höhna et al. (2022), we may find very distinct congruent scenarios, some of them not very plausible biologically, with difficulties to reach conclusions about the true diversification dynamics. Second, these sampling schemes seem to have difficulties sampling diversification histories close to the one used to define the congruence class. We certainly want to be able to sample different congruent trajectories, but recovery of the initial trajectories should be a prerequisite.

In this paper, we develop alternative sampling strategies of the congruence class, as well as a potential additional filtering step which consists in selecting the most parsimonious trajectories. We implement these approaches in CRABS. Finally, we illustrate their utility by applying them to a simulated scenario as well as mammals' diversification.

# 2 Methods

79

## 2.1 Sampling on the congruent curve

80

Given a congruence class uniquely characterized by the pulled diversification rate trajectory $r_p$ and the speciation rate at present $\lambda_0$, we aim to simultaneously sample speciation ($\lambda$) and extinction ($\mu$) rate trajectories within the class, instead of constraining

81
82
the functional space of one of the two trajectories. Following Höhna et al. (2022), we discretize time into small steps of size $\Delta t$ and use the linear approximation of $r_p$ mentioned in the Introduction. For each iteration $i$, we define the congruence curve as the set of $(\lambda_i, \mu_i)$ pairs such that the resulting trajectories remain within the congruence class. Rearranging the linear equation for $r_p$, the congruence curve has the following form (see Figure 1):

$$ \mu_i(\lambda_i) = \lambda_i - r_p(t_i) + \frac{1}{\lambda_i} \frac{\lambda_i - \lambda_{i-1}}{\Delta t} \quad (1) $$

At present $t_0$, $\lambda_0$ is given and we sample $\mu_0$ in a given distribution (see details below). Next, going backwards in time, we sample $\lambda_i$ and $\mu_i$ simultaneously along the congruence curve at each iteration $i$. This allows for a more flexible and balanced sampling within the congruence class.

Before sampling on the congruent curve, we compute two useful quantities: $\lambda_i^*$, the $\lambda_i$ value such as $\mu_i = \mu_{i-1}$, and $\lambda_i^{min}$, the minimum $\lambda_i$ value that maintains non-negative $\mu_i$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Sampling the congruence curve. Space of parameters $(\lambda_i, \mu_i)$ that can be sampled at a given iteration $i$. At each step, $\lambda_{i-1}$ and $r_p(t_i)$ define a specific congruence curve (in brown) along which $(\lambda_i, \mu_i)$ can be sampled so that the trajectory remains in the congruence class; the fact that it is below the 1:1 gray line in this illustration indicates a positive net diversification trend. We further impose the constraint of a minimal $\lambda_{min}$ value to avoid negative extinction rates (dashed area). Note that choosing $\lambda_i = \lambda_{i-1}$ implies a constant speciation rate while choosing $\lambda_i = \lambda_i^*$ implies a constant extinction rate; we will be sampling a position on the curve relative to these two values. The figure represents the situation when $\lambda_i^* > \lambda_{i-1}$, in which case $\lambda_i > \lambda_{i-1}$ and $\mu_i < \mu_{i-1}$, but we can also later have $\lambda_i^* < \lambda_{i-1}$, in which case $\lambda_i < \lambda_{i-1}$ and $\mu_i > \mu_{i-1}$.</figcaption></figure>
92 Assuming that $\forall i, \lambda_i > 0$:

$$ 
\begin{aligned}
\mu_i(\lambda_i) = \mu_{i-1} &\Leftrightarrow \lambda_i - r_p(t_i) - \frac{1}{\lambda_i} \frac{\lambda_i - \lambda_{i-1}}{\Delta t} = \mu_{i-1} \\
&\Leftrightarrow \lambda_i^2 + \left( \frac{1}{\Delta t} - r_p(t_i) - \mu_{i-1} \right) \lambda_i - \frac{\lambda_{i-1}}{\Delta t} = 0
\end{aligned}
$$

93 The only positive solution of this trinomial (given that the congruent curve is monotonic for $\lambda_i > 0$) is the following:

$$ \lambda_i^* = \frac{r_p(t_i) + \mu_{i-1} - \frac{1}{\Delta t} + \sqrt{(r_p(t_i) + \mu_{i-1} - \frac{1}{\Delta t})^2 + 4 \frac{\lambda_{i-1}}{\Delta t}}}{2} \quad (2) $$

94 From this equation, we can directly deduce $\lambda_i^{min}$, the minimum $\lambda_i$ value that keeps $\mu_i$ non-negative, by setting $\mu_{i-1}$ to 0:

$$ \lambda_i^{min} = \frac{r_p(t_i) - \frac{1}{\Delta t} + \sqrt{(r_p(t_i) - \frac{1}{\Delta t})^2 + 4 \frac{\lambda_{i-1}}{\Delta t}}}{2} \quad (3) $$

95 The iterative sampling along the congruence curve is illustrated in Figure 1, plotted with Geogebra (Hohenwarter, 2023).

## 2.2 Sampling approaches

97 At a given iteration $i$, on one extreme $\lambda$ is set constant ($\lambda_i = \lambda_{i-1}$), implying that all the information in the tree will be reflected in $\mu_i$. On the other extreme $\mu$ is set constant ($\lambda_i = \lambda_i^*$) and $\lambda_i$ conveys all the information. In order to obtain other sampling approaches, we select a position $p_i$ relative to $\lambda_{i-1}$ and $\lambda_i^*$. The various proposed approaches make different choices for selecting $p_i$. Below, "p-trajectory" refers to the trajectory of $p_i$ as we move through time steps $t_i$ from the present to the past.

### 2.2.1 Fixed position on the congruence curve: Constant p-trajectory

102 Our first approach consists in taking at each iteration $i$ the same relative position $p_0$ between $\lambda_{i-1}$ and $\lambda_i^*$ (i.e. $p_i$ is constant). Here is the complete procedure:

104 1. Choice of the position : $p_0 \sim Beta(\beta_1, \beta_2)$

105 2. Initialization : 
$$ \begin{cases} \lambda_0' = \lambda_0 \\ \mu_0' \sim \log \mathcal{N}(\log(\mu_0), \log \sigma_\mu) \end{cases} $$

106 3. For each iteration $i \ge 1$

107 (a) Compute $\lambda_i^{min}$ and $\lambda_i^*$

108 (b) Take $\lambda_i' = \lambda_{i-1} + p_0 \times (\lambda_i^* - \lambda_{i-1})$, i.e. at a relative position $p_0$ between $\lambda_{i-1}$ and $\lambda_i^*$

109 (c) If $\lambda_i' < \lambda_i^{min}$, set $(\lambda_i', \mu_i') = (\lambda_i^{min}, 0)$ to avoid negative $\mu$ values
(d) Else, compute $\mu'_{i} = \lambda'_{i} - r_{p}(t_{i}) + \frac{1}{\lambda'_{i}} \frac{\lambda'_{i} - \lambda'_{i-1}}{\Delta t}$ according to the congruence relationship

Note that $p_0 = 0$ implies a constant speciation rate throughout a clade's history (barring negative extinction rates), while $p_0 = 1$ implies a constant extinction rate. Our method explores intermediate trajectories between these two poles. As a consequence, even if the algorithm above only samples $\lambda$ before retrieving $\mu$ through the congruent curve, $\mu$ is still included in the procedure through the computation of $\lambda_i^*$. Note also that at each step this procedure forces either speciation rates to increase and extinction rates to decrease (if $\lambda_{i-1} < \lambda_i^*$), or the reverse (if $\lambda_{i-1} > \lambda_i^*$), which is in itself a constraint favoring parsimonious trajectories (as the data could typically also be explained by arbitrary large simultaneous increases in both rates).

The parameters of the Beta distribution $(\beta_1, \beta_2)$ can be set to $(1, 1)$ to sample uniformly between 0 and 1, or changed to different values in order to get more samples of specific trajectories (e.g. nearly constant lambda or mu, and all intermediates).

### 2.2.2 Autocorrelated position on the congruence curve: Markov Random Field (MRF) p-trajectory

This second approach is very similar to the previous one, with the distinction that $p_i$ is not constant at each iteration $i$, but instead follows an autocorrelated process starting at $p_0$. Following Magee et al. (2020), we use either a Gaussian Markov random field (GMRF), a model in which $p_i|(p_{i-1}, \gamma) \sim N(p_{i-1}, \gamma^2)$, or a Horseshoe Markov random field (HSMRF), a model in which $p_i|(p_{i-1}, \zeta) \sim Horseshoe(p_{i-1}, \gamma^2)$. $\gamma$ is a global scale parameter that controls the smoothness of the overall field (shrinking all jumps towards zero when $\gamma$ tends towards zero). In the case of the HSMRF, the fat-tailed nature of the Horseshoe distribution allows for rare occurrences of larger jumps in the trajectory.

1. Drawing the p-trajectory

(a) **Global scale parameter** controlling the degree of regularization (i.e. steering towards "simpler" trajectories) : $\gamma \sim \text{halfCauchy}(0, \zeta)$

This parameter is itself controlled hierarchically by the hyperparameter $\zeta$, which can be obtained with the CRABS function `get.gmrf.global.scale` such that the expected number of 2-fold $p$ changes given the number of time steps equals $\log(2)$.

It can be further adjusted manually with a scaling factor $\sigma_{\text{MRF}}$.

(b) Trajectory: starting value $p_0 \sim \text{Beta}(\beta_1, \beta_2)$, then

$$ \begin{cases} \text{GMRF: } p_{i \ge 1} \sim N(p_{i-1}, sd = \sigma_{\text{MRF}} \cdot \gamma) \\ \text{HSMRF: } p_{i \ge 1} \sim N(p_{i-1}, sd = \sigma_{\text{MRF}} \cdot \sigma_i \cdot \gamma) \text{ with } \sigma_i \sim \text{halfCauchy}(0, 1) \end{cases} $$

$\sigma_i$ is a **local shrinkage parameter** allowing for rare larger jumps.

(c) Centering: by default, such autocorrelated trajectories will see their variance increase over the time steps, resulting in more variable dynamics in the past than in the present. To cancel this effect, the p-trajectory is centered around the $p_i$ value at a randomly sampled time step.
(d) Rescaling: here the p-trajectories are not necessarily bounded between 0 and 1, so it is possible for $\lambda_i$ and $\mu_i$ to both increase ($p_{max} > 1$) or decrease ($p_{min} < 0$) simultaneously. We however rescale the p-trajectories whenever they leave predefined bounds. We used $p_{min} = -0.05$ and $p_{max} = 1.05$, which worked well in practice, and we implemented them as default values.

2. Initialization: $$ \begin{cases} \lambda_0' = \lambda_0 \\ \mu_0' \sim \log N(\log(\mu_0), \log \sigma_\mu) \end{cases} $$
3. For each iteration $i \ge 1$: same as before with $p_i$ instead of $p_0$

### 2.2.3 MRF p-trajectory with a global regularity filter

It is possible to construct congruent diversification rate trajectories with arbitrarily large and sudden rate changes, however these diversification histories may be less plausible biologically than smoother ones. If one follows the parsimony principle, simple trajectories that explain the observed data as well as more complex ones should be preferred Morlon et al. (2022). This a priori intuition can be incorporated by generating many congruence diversification trajectories, for example following the MRF p sampling approach, and selecting those that do not exhibit such sudden sharp changes, for example using a regularity criterion that can be chosen freely. We tested here three regularity criteria:

* $\mathcal{L}_1$ norm on all the step-wise variations in both speciation ($\lambda_i$) and extinction rates ($\mu_i$), for $i \in [0; N]$, where N is the number of time steps:

$$ \mathcal{L}_1 = \sum_{i=1}^{N} [abs(\lambda_i - \lambda_{i-1}) + abs(\mu_i - \mu_{i-1})] $$

* $\mathcal{L}_2$ norm on step-wise variations:

$$ \mathcal{L}_2 = \sqrt{\sum_{i=1}^{N} [(\lambda_i - \lambda_{i-1})^2 + (\mu_i - \mu_{i-1})^2]} $$

* second-order $\mathcal{L}_1$ norm on derivative shifts in both speciation and extinction rates:

$$ \mathcal{L}_1' = \sum_{i=2}^{N} [abs(\lambda_i - 2\lambda_{i-1} + \lambda_{i-2}) + abs(\mu_i - 2\mu_{i-1} + \mu_{i-2})] $$

### 2.2.4 Implementation in CRABS

We implemented our new sampling methods in R as extended functions of the CRABS package (R Core Team, 2018; Höhna et al., 2022). As in CRABS, two functions are called. First, the function `sample.basicConfig.joint` (modified from `sample.basicConfig`), which in its default setting returns a sampler for HSMRF p-trajectories, and can return a sampler for constant p-trajectories by choosing a zero scaling factor (`mrf.sd.scale = 0`). We chose the HSMRF strategy as a default, since the GMRF is less flexible. Second, the function `sample.congruence.class` returns the samples; we simply added the argument `sample.joint.rates` to this function, and the user should specify to use the option `rate.type = "joint"`.
We also implemented the function `full.plot.regularity threshholds` that ranks the sampled trajectories according to their regularity values, selects the x% most regular, and displays the corresponding trajectories. The arguments `filtering_fractions` (by default 1%, 5%, 20% and 90%), `penalty` (among "L1", "L2" and "L1_derivative" for the derivative shifts regularity criterion) and `rates` (vector including "lambda" for speciation, "mu" for extinction, "delta" for net-diversification, "epsilon" for turnover) give control over the specific filtering being performed.

## 2.3 Illustrations

### 2.3.1 Simulated scenario

In order to illustrate the utility of our sampling approach in comparison with previous approaches, we analyse a mock reference scenario with a declining speciation rate ($\lambda(t) = 0.1 + \exp(\frac{t-70}{10})$) and a mostly constant extinction rate with a sudden spike ($\mu(t) = 0.05 + f_{\beta}(\frac{t}{30}, 100, 100)$) / 20 with $f_{\beta}$ the probability density function of the Beta distribution). This model corresponds to an often-observed scenario of early-burst diversification (Moen & Morlon, 2014) combined with a mass extinction event.

We first apply the previous CRABS functions (Höhna et al., 2022), sampling 20 congruent HSMRF trajectories of either speciation or extinction rates using the default options, and summarize relative rate changes over time with the recommended threshold of 0.02 for significant increases and decreases.

We then carry out our two new procedures on the same scenario, sampling 20 constant p-trajectories and 20 HSMRF p-trajectories, each corresponding to a joint sampling of speciation and extinction rates.

### 2.3.2 Mammalian diversification

We analyse the diversification dynamics of mammals, using the dated mammalian phylogeny constructed by Álvarez-Carretero et al. (2022), which integrates information from 72 mammal genomes. More specifically, we used the subset of this timetree that excludes subspecies (634 among the 4705 initial lineages), as used in Quintero et al. (2022).

We conducted an initial diversification analysis in RevBayes (Höhna et al., 2016) assuming auto-correlated piecewise-constant speciation and extinction rates trajectories described by a HSMRF prior distribution (see Magee et al. (2020) and the corresponding RevBayes tutorial<fnref n="1" />). We used 30 time intervals whose limits were determined during the inference process. We applied clade-specific sampling fractions, using the number of missing species for 16 monophyletic subclades provided by Quintero et al. (2022). We visualized the results using the RevGadgets R library (Tribble et al., 2022; R Core Team, 2018).

Finally, we generated 5000 joint HSMRF p-trajectories on the congruent curve, and filtered them based on their regularity values, using each of the three available norms. For illustration purposes, we sampled four sets of 50 trajectories within those that fell below the 1st, 5th, 20th, and 90th centiles of regularity values. The number of initial samples must be high enough to ensure that the 1st percentile subsampling represents enough trajectories. We therefore recommend generating at least 100 times the number of retained samples for each selection threshold (50 in our case).

<footnote n="1">¹https://revbayes.github.io/tutorials/divrate/ebd.html</footnote>
187

# 3 Results

188

## 3.1 Simulated Scenario

189

Applied to the mock scenario with a decline in speciation rates and a spike in extinction, the methods of Höhna et al. (2022), which sample either speciation or extinction rate trajectories, lead to two possibilities depending on whether one or the other is sampled (Fig. 2). Sampling the extinction rate trajectory allows recovering the decline in speciation rate, but also infers a systematic short increase in the speciation rate at the time of the (unaccounted for) extinction peak (Fig 2 a,c). Sampling the speciation rate trajectory instead allows recovering the extinction peak, but also infers a spurious initial increase in extinction rate reflecting the (unaccounted for) decline in speciation (Fig 2 b,d). In addition, all inferred extinction rates at the origin of the clade are negative, illustrating the difficulty of sampling the congruence class through speciation trajectories.

190

191

192

193

194

195

196

In contrast, the sampling approaches we developed were able to sample a greater diversity of congruent scenarios, including around the initial model (Fig. 3). With the constant *p*-trajectory sampling approach (Fig. 3a), extinction rates ranged from completely flat to the large initial spike, and speciation rates showed a steady decline, with sometimes a small decrease at the time of the sharp decline in the net diversification rate, corresponding to the cases when the spike is not fully reproduced in the extinction curve. HSMRF *p*-trajectories (Fig. 3b) result in similar sampled trajectories with more variability, including congruent scenarios with higher speciation and extinction rates in the past.

197

198

199

200

201
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

**Figure 2: Illustration of current limitations in the exploration of congruence classes.** (a,b) Rate trajectories over time for the reference scenario (black line) and 20 congruent trajectories (dashed lines) obtained with previous CRABS functions (Höhna et al., 2022) by sampling HSMRF extinction (a) or speciation (b) rate trajectories. (c,d) Top: relative rate changes over time, with red dashed lines demarcating the 0.02 significance threshold. Bottom: summary across 10 congruent models with increases in green and decreases in purple; the first row at the top corresponds to the reference scenario. Note that in both cases, no sampled congruent trajectory shows patterns close to the reference scenario used to construct the congruence class, and that the macroevolutionary signal is reflected either entirely on speciation rate, or entirely on extinction rate.
<figure><img src="image_7.png" /><figcaption>(b) *HSMRF p trajectories.*</figcaption></figure>

<figure><img src="image_8.png" /></figure>

**Figure 3: Two approaches to jointly sample speciation rate and extinction rate trajectories.** Each blue *p*-trajectory on the left panel corresponds to rate trajectories on the right panel. The black line represents the mock diversification scenario used in Figure 2, other trajectories (dashed lines) are sampled from its congruence class. Their trends are summarized in Sup. Fig. A.3. (a) Evolutionary rates are sampled through time always at the same position between between $\lambda_{i-1}$ ($\lambda$ as constant as possible without negative extinction) and $\lambda_i^*$ (constant $\mu$). (b) Evolutionary rates are sampled through time similarly, except that the position relative to $\lambda_{i-1}$ and $\lambda_i^*$ follows a HSMRF process.

## 3.2 Mammalian diversification histories

The diversification analyses on the mammalian phylogeny assuming auto-correlated piecewise-constant speciation and extinction rates trajectories in RevBayes suggest a 10-fold jump in speciation rates (from 0.01 to $0.1 \text{ Msy}^{-1}$) some 80-85 million years ago and a constant, low extinction rate (Figure 4). Congruent mammalian diversification histories that are selected under the most stringent regularity constraints (the 1% most regular trajectories) match closely the initial trajectories (Figure 4).

Some of the more complex trajectories paint a different picture, where the 80-85 Mya signal is driven both by fluctuations in speciation rates and by increased extinction rates before and during this period. These congruent trajectories are also characterised by markedly lower diversification rates during the transition period. Similar patterns are found with other regularity penalties (see Supplementary Figures B.1 and B.2).
<figure><img src="image_9.png" /></figure>

Figure 4: Congruent mammalian diversification trajectories for increasing regularity thresholds. The median rate trajectory inferred in RevBayes is displayed in black, together with congruent speciation (top), extinction (middle) and diversification (bottom) rates in dashed lines. Joint HSRMF *p*-trajectories were generated on the congruent curve, and selected based on their regularity. The most stringent regularity criterion is represented on the left, with 50 trajectories sampled below the first centile of regularity values. Similar plots are provided using the 5th, 20th, and 90th centiles, from left to right.

# 4 Discussion

211

212

213

214

215

216

217

218

219

We introduced a novel approach for sampling alternative diversification scenarios in a given congruence class. By jointly sampling congruent speciation and extinction rate trajectories, our approach provides a more thorough representation of congruent classes than previous methods which sampled entirely one of the two trajectories before deducing the other. Although a congruent class contains very distinct diversification scenarios, which are indistinguishable based on their likelihood to have produced a given extant phylogeny, these scenarios are not equally regular. We implemented the possibility to select congruent scenarios based on several regularity criteria. As we illustrated in our empirical analysis of mammals, this filtering step can shrink the congruence class to a unique diversification scenario. Our developments therefore provide the possibility to explore the full range of equally likely diversification scenarios while keeping biological interpretability.

## 4.1 Embracing the uncertainty

220

221

222

223

224

225

In the face of the non-identifiability of birth-death diversification models with unconstrained time-variable rates, it has been proposed to restrict inferences to identifiable reparametrizations of the model Louca & Pennell (2020). Unfortunately, such reparametrizations come at the cost of limited biological interpretability (Helmstetter et al., 2022). Another approach is to work with hypotheses-driven models (Morlon et al., 2022). In this case the models often are identifiable, but this might give overconfidence in misleading results, as the conclusions may be highly dependent on the hypotheses being made. These two
approaches have the common disadvantage to force us to restrict our analysis to identifiable patterns, rather than to formulate
the models we believe best reflect the processes at play regardless of whether they are identifiable or not (Gustafson et al.,
2005; Gustafson, 2015). Yet another approach is to rely on a Bayesian framework, in which inference of posterior parameter
distributions is always theoretically identifiable. However, when the underlying model is non-identifiable, the final outcomes of
Bayesian inference are highly dependent on the choice of priors, reflecting our current beliefs. In such situations, it can be much
beneficial to embrace the uncertainty and complement the analyses with an exploration of congruent trajectories, which represent
a less subjective probing of the potential processes at work. Such exploration can help assess the results' sensitivity to the a
*priory* hypotheses and/or the prior choice: congruent trajectories that deviate significantly from the originally inferred posterior
trajectories but appear biologically realistic suggest that other hypotheses or priors may warrant further scrutiny.

Höhna et al. (2022) paved the way to the exploration of congruence classes, yet their sampling approach is not as exhaustive as it
may seem. Sampling one rate trajectory arbitrarily as a first step before computing the other trajectory to obtain a congruent
scenario in a second step forces signal in the data to be captured in the second trajectory. One risk, given the difficulty of
obtaining trajectories with positive extinction rates when sampling speciation rate trajectories, is to sample mostly extinction rate
trajectories (see e.g. Kopperud et al. (2023)), and therefore to force signal in the data to be entirely captured in the speciation rate
trajectory. This provides a biased exploration of congruence classes that may give an impression of robustness that could fall
apart when exploring the class more exhaustively. Sampling congruent speciation and extinction rate trajectories jointly, as we
do here, provides an efficient approach for exploring the congruence class while allowing signal in the data to be captured by
both speciation and extinction rate curves, therefore truly embracing the uncertainty.

## 4.2 Shrinking the uncertainty

Although embracing the uncertainty by exploring congruence classes provides a fair representation of the diversity of scenarios
that cannot be distinguished from their likelihood, it can be useful to shrink this uncertainty using other criteria, such as regularity
criteria. Facing equally likely hypotheses is not uncommon in the process of scientific inquiry, but this has never been an
insurmountable barrier in the face of Occam's razor blade. Indeed, it is always possible to create more complex models to
explain any set of observations, but hypotheses with equal likelihoods need not be equally plausible. Among a set of equally
likely hypotheses, it is common practice to assign higher prior probabilities to the most parsimonious ones (Morlon et al., 2022).
We have implemented several regularity criteria to perform this selection. Interestingly, in our empirical illustration on mammals,
the most regular congruent diversification scenarios match closely the initial scenario obtained with a Bayesian implementation
of an auto-correlated piecewise-constant diversification model, suggesting that currently-used diversification models already
perform well in selecting the most regular trajectories within a congruence class.

Visualizing congruent rates under varying stringency thresholds (as depicted in Fig. 4) eliminates the need to make an arbitrary
choice for the threshold and transparently displays the distribution of possible trajectories, leaving the interpretation to depend
on the reader's willingness to favor the "simplest" histories that can explain the data.
**4.3 Perspectives**

Congruence classes have been characterized beyond classical birth-death models, for example in the context of the fossilized birth-death process. This process contains an additional rate for incorporating samples over time into the phylogeny, interpreted as fossils in macroevolutionary studies and genetic samples in epidemiological analyses (Andréoletti et al., 2022; MacPherson et al., 2022). Although this model has been demonstrated to be non-identifiable in the (macroevolutionarily unrealistic) case where lineages go extinct upon sampling (Louca et al., 2021), a more general demonstration is lacking and a method for sampling within this congruence class has not been established. The methodology presented in this study could potentially be extended for the joint sampling of the three rate parameters.

**5 Conclusion**

Embracing the uncertainty inherent to phylogenetic diversification analyses by exploring congruence classes can help us avoid missing realistic diversification scenarios that cannot be ruled out by the observed data. As we nevertheless strive to narrow down this uncertainty, we contend that restricting ourselves to simpler identifiable quantities that are hard to interpret biologically is not the most promising strategy. Instead, introducing additional independent data types and prior knowledge while navigating between embracing and shrinking the uncertainty, can provide a reliable and nuanced portrayal of past diversification histories, ultimately enhancing our comprehension of the complex evolutionary processes that shape ecological systems.

**6 Acknowledgements**

We thank Bjørn T. Kopperud, Sebastian Höhna, Andrew F. Magee for their constructive feedback and reviewing the code, and
Lucas Rey and Nathanaël Boutillon for discussions on the sampling method.

**7 Competing interests**

The authors report no conflict of interest.

**8 Availability of data and materials**

Our functions have been integrated into the CRABS R package, available on Github (https://github.com/afmagee/CRABS).
The code and data corresponding to the analyses presented in this article are available at https://doi.org/10.5281/zenodo.
8171427.
# 9 Funding information

J.A. is supported by a doctoral fellowship from the École Normale Supérieure de Paris. HM acknowledges support from the

European Research Council (grant CoG-PANDA).

# 10 Authors' contributions

Both authors contributed to the design of the project, the discussion of the research, and the writing and editing of the article.

Jérémy Andréoletti wrote the code, performed the analysis and wrote the first draft.
## References

288

289

290

291

Helmstetter, A., Glemin, S., Käfer, J., Zenil-Ferguson, R., Sauquet, H., Boer, H., Dagallier, L., Mazet, N., Reboud, E., Couvreur, T. & Condamine, F. Pulled Diversification Rates, Lineages-Through-Time Plots, and Modern Macroevolutionary Modeling. *Systematic Biology*. **71**, 758-773 (2022,5), https://doi.org/10.1093/sysbio/syab083

292

293

Louca, S. & Pennell, M. Extant timetrees are consistent with a myriad of diversification histories. *Nature*. **580**, 502-505 (2020,4), https://www.nature.com/articles/s41586-020-2176-1, Number: 7804 Publisher: Nature Publishing Group

294

295

296

Höhna, S., Kopperud, B. & Magee, A. CRABS: Congruent rate analyses in birth–death scenarios. *Methods In Ecology And Evolution*. n/a (2022), https://onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.13997, eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/2041-210X.13997

297

298

299

Quintero, I., Lartillot, N. & Morlon, H. The birth-death diffusion leading to present-day Mammal diversity. (bioRxiv,2022,8), https://www.biorxiv.org/content/10.1101/2022.08.09.503355v1, Pages: 2022.08.09.503355 Section: New Results

300

301

302

Álvarez-Carretero, S., Tamuri, A., Battini, M., Nascimento, F., Carlisle, E., Asher, R., Yang, Z., Donoghue, P. & Reis, M. A species-level timeline of mammal evolution integrating phylogenomic data. *Nature*. **602**, 263-267 (2022,2), https://www.nature.com/articles/s41586-021-04341-1, Number: 7896 Publisher: Nature Publishing Group

303

304

305

Höhna, S., Landis, M., Heath, T., Boussau, B., Lartillot, N., Moore, B., Huelsenbeck, J. & Ronquist, F. RevBayes: Bayesian Phylogenetic Inference Using Graphical Models and an Interactive Model-Specification Language. *Systematic Biology*. **65**, 726-736 (2016,7)

306

307

308

Magee, A., Höhna, S., Vasylyeva, T., Leaché, A. & Minin, V. Locally adaptive Bayesian birth-death model successfully detects slow and rapid rate shifts. *PLOS Computational Biology*. **16**, e1007999 (2020,10), https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007999, Publisher: Public Library of Science

309

310

311

Tribble, C., Freyman, W., Landis, M., Lim, J., Barido-Sottani, J., Kopperud, B., Höhna, S. & May, M. RevGadgets: An R package for visualizing Bayesian phylogenetic analyses from RevBayes. *Methods In Ecology And Evolution*. **13**, 314-323 (2022,2), https://onlinelibrary.wiley.com/doi/10.1111/2041-210X.13750

312

313

R Core Team R: A language and environment for statistical computing.. (R Foundation for Statistical Computing, Vienna, Austria.,2018,1), https://www.R-project.org/

314

315

Louca, S. & Doebeli, M. Efficient comparative phylogenetics on large trees. *Bioinformatics*. **34**, 1053-1055 (2018,3), https://academic.oup.com/bioinformatics/article/34/6/1053/4582279

316

317

318

Louca, S., McLaughlin, A., MacPherson, A., Joy, J. & Pennell, M. Fundamental Identifiability Limits in Molecular Epidemiology. *Molecular Biology And Evolution*. **38**, 4010-4024 (2021,8), https://academic.oup.com/mbe/article/38/9/4010/6278301
319

Lambert, A. & Stadler, T. Birth–death models and coalescent point processes: The shape and probability of reconstructed phylogenies. *Theoretical Population Biology*. **90**, pp. 113-128 (2013,12), https://www.sciencedirect.com/science/article/pii/S0040580913001056

320

321

322

Stadler, T. On incomplete sampling under birth–death models and connections to the sampling-based coalescent. *Journal Of Theoretical Biology*. **261**, 58-66 (2009,11), https://www.sciencedirect.com/science/article/pii/S0022519309003300

323

324

325

Morlon, H., Robin, S. & Hartig, F. Studying speciation and extinction dynamics from phylogenies: addressing identifiability issues. *Trends In Ecology & Evolution*. **37**, 497-506 (2022,6), https://www.sciencedirect.com/science/article/pii/S0169534722000271

326

327

328

Stadler, T. How Can We Improve Accuracy of Macroevolutionary Rate Estimates?. *Systematic Biology*. **62**, 321-329 (2013,3), https://academic.oup.com/sysbio/article/62/2/321/1667317

329

330

Legried, B. & Terhorst, J. A class of identifiable phylogenetic birth-death models. *Proceedings Of The National Academy Of Sciences Of The United States Of America*. **119**, e2119513119 (2022,8), https://www.pnas.org/doi/10.1073/pnas.2119513119

331

332

333

Legried, B. & Terhorst, J. Identifiability and inference of phylogenetic birth-death models. (bioRxiv,2022,10), https://www.biorxiv.org/content/10.1101/2022.08.26.505438v2, Pages: 2022.08.26.505438 Section: New Results

334

335

Raup, D. Cohort analysis of generic survivorship. *Paleobiology*. **4**, 1-15 (1978), https://www.cambridge.org/core/journals/paleobiology/article/abs/cohort-analysis-of-generic-survivorship/94FFA08A22C590776CECD137D2195602, Publisher: Cambridge University Press

336

337

338

Foote, M. Temporal variation in extinction risk and temporal scaling of extinction metrics. *Paleobiology*. **20**, 424-444 (1994), https://www.cambridge.org/core/journals/paleobiology/article/abs/temporal-variation-in-extinction-risk-and-temporal-scaling-of-extinction-metrics/9A70EDB2AF955F38B133AF6E9F3F8388, Publisher: Cambridge University Press

339

340

341

342

Foote, M. Origination and extinction components of taxonomic diversity: General problems. *Paleobiology*. **26**, 74-102 (2000), https://www.cambridge.org/core/journals/paleobiology/article/abs/origination-and-extinction-components-of-taxonomic-diversity-general-problems/7FE12B43106FC20C9AFC4031F72A56F9, Publisher: Cambridge University Press

343

344

345

346

Jablonski, D. & Chaloner, W. Extinctions in the Fossil Record [and Discussion]. *Philosophical Transactions: Biological Sciences*. **344**, 11-17 (1994), https://www.jstor.org/stable/56148, Publisher: The Royal Society

347

348

Alroy, J. Accurate and precise estimates of origination and extinction rates. *Paleobiology*. **40**, 374-397 (2014), https://www.cambridge.org/core/journals/paleobiology/article/abs/accurate-and-precise-estimates-of-origination-and-extinction-rates/AB73EA51B6748D3843199A8A0D6B1603, Publisher: Cambridge University Press

349

350

351
352

Silvestro, D., Schnitzler, J., Liow, L., Antonelli, A. & Salamin, N. Bayesian Estimation of Speciation and Extinction from Incomplete Fossil Occurrence Data. *Systematic Biology*. **63**, 349-367 (2014,5), https://doi.org/10.1093/sysbio/syu006

353

354

355

Nee, S., May, R. & Harvey, P. The reconstructed evolutionary process. *Philosophical Transactions Of The Royal Society Of London. Series B: Biological Sciences*. **344**, 305-311 (1994,5), https://royalsocietypublishing.org/doi/abs/10.1098/rstb.1994.0068, Publisher: Royal Society

356

357

358

Morlon, H., Parsons, T. & Plotkin, J. Reconciling molecular phylogenies with the fossil record. *Proceedings Of The National Academy Of Sciences*. **108**, 16327-16332 (2011,9), https://www.pnas.org/doi/full/10.1073/pnas.1102543108, Publisher: Proceedings of the National Academy of Sciences

359

360

361

Stadler, T. Sampling-through-time in birth–death trees. *Journal Of Theoretical Biology*. **267**, 396-404 (2010,12), https://www.sciencedirect.com/science/article/pii/S0022519310004765

362

363

Alfaro, M., Santini, F., Brock, C., Alamillo, H., Dornburg, A., Rabosky, D., Carnevale, G. & Harmon, L. Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates. *Proceedings Of The National Academy Of Sciences*. **106**, 13410-13414 (2009,8), https://www.pnas.org/doi/full/10.1073/pnas.0811087106, Publisher: Proceedings of the National Academy of Sciences

364

365

366

367

Moen, D. & Morlon, H. Why does diversification slow down?. *Trends In Ecology & Evolution*. **29**, 190-197 (2014,4), https://www.sciencedirect.com/science/article/pii/S016953471400024X

368

369

Hohenwarter, M. Geogebra 6. (2023,1), http://www.geogebra.org

370

Kopperud, B., Magee, A. & Höhna, S. Rapidly changing speciation and extinction rates can be inferred in spite of nonidentifiability. *Proceedings Of The National Academy Of Sciences*. **120**, e2208851120 (2023,7), https://www.pnas.org/doi/abs/10.1073/pnas.2208851120, Publisher: Proceedings of the National Academy of Sciences

371

372

373

Lan, S., Palacios, J., Karcher, M., Minin, V. & Shahbaba, B. An efficient Bayesian inference framework for coalescent-based nonparametric phylodynamics. *Bioinformatics*. **31**, 3282-3289 (2015,10), https://doi.org/10.1093/bioinformatics/btv378

374

375

376

Gustafson, P., Gelfand, A., Sahu, S., Johnson, W., Hanson, T., Joseph, L. & Lee, J. On Model Expansion, Model Contraction, Identifiability and Prior Information: Two Illustrative Scenarios Involving Mismeasured Variables [with Comments and Rejoinder]. *Statistical Science*. **20**, 111-140 (2005), https://www.jstor.org/stable/20061166, Publisher: InstituteofMathematicalStatistics

377

378

379

380

Andréoletti, J., Zwaans, A., Warnock, R., Aguirre-Fernández, G., Barido-Sottani, J., Gupta, A., Stadler, T. & Manceau, M. The Occurrence Birth–Death Process for Combined-Evidence Analysis in Macroevolution and Epidemiology. *Systematic Biology*. **71**, 1440-1452 (2022,11), https://doi.org/10.1093/sysbio/syac037

381

382

383

MacPherson, A., Louca, S., McLaughlin, A., Joy, J. & Pennell, M. Unifying Phylogenetic Birth–Death Models in Epidemiology and Macroevolution. *Systematic Biology*. **71**, 172-189 (2022,1), https://doi.org/10.1093/sysbio/syab049

384
385 Stadler, T. Mammalian phylogeny reveals recent diversification rate shifts. *Proceedings Of The National Academy Of Sciences*.
386 **108**, 6187-6192 (2011,4), https://www.pnas.org/doi/10.1073/pnas.1016876108, Publisher: Proceedings of the
387 National Academy of Sciences

388 Rabosky, D. Automatic Detection of Key Innovations, Rate Shifts, and Diversity-Dependence on Phylogenetic Trees. *PLOS ONE*. **9**, e89543 (2014,2), https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0089543,
389 Publisher: Public Library of Science
390

391 Barido-Sottani, J., Vaughan, T. & Stadler, T. A Multitype Birth–Death Model for Bayesian Inference of Lineage-Specific Birth
392 and Death Rates. *Systematic Biology*. **69**, 973-986 (2020,9), https://doi.org/10.1093/sysbio/syaa016

393 Maliet, O., Hartig, F. & Morlon, H. A model with many small shifts for estimating species-specific diversification rates. *Nature Ecology & Evolution*. **3**, 1086-1092 (2019,6), https://www.nature.com/articles/s41559-019-0908-0

395 Kubo, T. & Iwasa, Y. Inferring the Rates of Branching and Extinction from Molecular Phylogenies. *Evolution*. **49**, 694-704
396 (1995), https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1558-5646.1995.tb02306.x, eprint: https://
397 onlinelibrary.wiley.com/doi/pdf/10.1111/j.1558-5646.1995.tb02306.x

398 Gustafson, P. Bayesian Inference for Partially Identified Models: Exploring the Limits of Limited Data. (CRC Press,2015,4),
399 Google-Books-ID: p4y9BwAAQBAJ
# Appendix A Simulated scenario

## A.1 Exponential decrease of speciation rates or exponential increase of extinction rates

<figure><img src="image_10.png" /><figcaption>(b) Sampling speciation rate as exponential decreasing trajectories.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Summary of speciation trends.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) Summary of extinction trends.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

Figure A.1: Illustration of diversification scenarios recovered with current congruence class sampling methods. Same as Figure 2, with speciation or extinction rates sampled as exponential trajectories.

## A.2 GMRF trajectories of speciation and extinction rates
<figure><img src="image_14.png" /><figcaption>(c) Summary of speciation trends.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) Summary of extinction trends.</figcaption></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

**Figure A.2: Illustration of diversification scenarios recovered with current congruence class sampling methods.** Same as Figure 2, with speciation or extinction rates sampled as GMRF processes.
## A.3 Joint sampling with *p*-trajectories

<figure><img src="image_18.png" /><figcaption>Figure A.3: Summaries of the congruent diversification scenarios shown in Figure 3. The left panel summarizes speciation trends and the right panel extinction trends, both with a 2% threshold for significant increases or decreases.</figcaption></figure>

## Appendix B Mammalian diversification histories
<figure><img src="image_19.png" /><figcaption>**Figure B.1:** Congruent mammalian diversification trajectories for increasing regularity thresholds. 5000 joint HSRMF trajectories were sampled on the congruent curve, and then for each column 50 of them were selected among the most regular one (i.e. below the 1st, 5th, 20th, and 90th centiles, respectively). The chosen regularity criterion is a step-wise $L^2$ norm on all the step-wise variations in both speciation and extinction rates.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure B.2:** Congruent mammalian diversification trajectories for increasing regularity thresholds. 5000 joint HSRMF trajectories were sampled on the congruent curve, and then for each column 50 of them were selected among the most regular one (i.e. below the 1st, 5th, 20th, and 90th centiles, respectively). The chosen regularity criterion is a step-wise $L^1$ norm on all derivative shifts in both speciation and extinction rates.</figcaption></figure>