To an extent, this makes me think back to the problem highlighted by Louca and Pennell (2020) of there being, for any given phylogeny of extant species only, large "congruence classes" of disparate diversification scenarios that could produce the observed trees. Could it be that this is still a problem that is faced by neural networks, as well as the more classical likelihood methods? Perhaps GNNs and the architectures you outline here could excel at inferring more localized (e.g. branch-specific, or tree-heterogeneous) diversification dynamics, rather than estimating tree-wide parameters under a global model?

I suspect/hope you did these ancestral state reconstructions on time-calibrated phylogenies? If so, I would just clarify this here given these methods are only described in the supplement, and this methodological detail will have strong impacts on their outcome, as it does not make sense to conduct ASR on non-time-calibrated trees.

While I understand that it's a logical extension to represent the tree as a vector/matrix format as you've done here to obtain summary statistics, it should on principle be feasilble to do something similar using graph neural networks, which can learn tree/graph representations directly from the original phylogeny/graph. Have you considered taking such an approach?

What do the different colors correspond to here (the figure is missing a legend)? You state in the text that the black line is the expectation, only describe one more - the inference limit, though it's unclear what color this is. I suspect this is the single line in red, whereas the estimated approximate posterior distribution from cxt is in blue?

I think this is a fantastic contribution, and it makes absolute sense that accounting for evolutionary relationships using the attention mechanism would afford such a performance boost!

I do wonder though, taking this idea one step forward, how much might be lost by using a fixed species tree topology to bias clade-level attention, as opposed to using per-window trees to bias per-sequence attention. Obviously the latter would come with significant additional computational cost, but I can't help but suspect that meaningful signal is being lost by not accounting for genomic heterogeneity of evolutionary relationships.

But see Leroy et al., 2025 (https://doi.org/10.1101/2025.08.14.670341) that addresses this issue using an improved pooling operator. However, as they discuss in their preprint, the performance they achieve (exceeding that of the MLE) still likely does not represent a ceiling on their performance here, as the architecture is quite simple. Use of more sophisticated graph-based architectures including graph transformers (which combat oversmoothing and can more readily account for both local and global patterns) will likely increase this performance further.

Related to the pooling operator, I think large gains may come from the use of 1) edge weights in your GCN layers so that not all neighbors are treated equally by the message passing mechanism, and 2) alternative MPNN layer types, including use of the graph attention mechanism (i.e. GAT) or graph transformers, which use the attention mechanism to learn which neighbors are more "important." I suspect that even with simple mean-pooling, these alternative layer types will be much more performant and generalizable (e.g. from CRBD to BiSSE). In effect the GCN layers (particularly without using edge weights) is more akin to the CRBD in that it assumes uniform, homogeneous contribution by all neighbors to feature updates.

Again, I think this is a case for exploring the use of more general pooling operators, such as EdgePooling, etc, that might capture the relevant signal, but without imposing such rigid inductive biases on the architecture that could prove harmful to more general application.

Have you considered other pooling layers other than mean-pooling? For instance existing pooling operators like EdgePooling might confer similar benefits by retaining the temporal signal through iterative edge contraction. https://pytorch-geometric.readthedocs.io/en/latest/generated/torch_geometric.nn.pool.EdgePooling.html#torch_geometric.nn.pool.EdgePooling

Related to the above (regarding the number of GNN layers used), have you considered/looked into using either gate residual connections, or using something like jumping knowledge (https://arxiv.org/abs/1806.03536) between layers to mitigate oversmoothing in deeper networks which GNNs tend to be prone to? This could be another simple modification with outsized benefits.

I presume you're using the "base" GCN layer here? If so I'd be clear about this, since unlike the CNN and MLP, there is a massive diversity of MPNN layer types - its worth being very explicit about what you're using, since their choice could have massive impacts on the performance of these architectures!

Is there a reason you chose not to include edge weights/attributes corresponding to each respective branch length? I suspect this would be quite useful/informative for providing additional context to the message passing layers and could boost performance further. Otherwise, each neighbor is assumed to contribute equally to feature updates, which is not something we innately suspect to be true.

True, but only necessarily when the prediction task is graph-level. This is not necessary for many other common GNN prediction tasks, including node or link prediction. This is a bit of a nitpick, but I think useful to distinguish since this architecture is so new to the field!

I know these details are present in the supplement, but here and for the other NNs, it could be useful to specify the number of layers used - particularly for the GNN, where this relates to the effective "diameter" of visibility or neighborhood size (i.e. how many hops away the message passing mechanism aggregates neighbor information).

It seems worth communicating here that in 5 (6?) out of 8 total antibiotics tested, the 2D-VAE outperformed the 2D-RHVAE for out-of-sample prediction. That said, it is noteworthy that, at least for one antibiotic (NQO), the magnitude of improvement over the VAE by the RHVAE is the largest observed.

I should emphasize too - I don't think this is necessarily a problem for the RHVAE - it clearly has significant benefits outside of the capacity to generalize! If anything, I'd say it's reasonable to hypothesize that perhaps it doesn't generalize as well because it does a better job of learning where it can and can not make accurate predictions, leading to sharper boundaries in prediction accuracy.

I'm very intrigued by this idea. At the very least, this seems like an interesting and effective way to drive VAEs towards learning more biologically meaningful, and increasingly expressive prior distributions rather than using a simple normal. On the other hand, this suggests the promise of a latent space that may even be interpretable.

I'm wondering if you've looked into this - whether these meaningful distances in the high-dimensional latent space are, or eventually could be made to be interpretable such that the contributions of different phenotypes to each latent dimension is quantifiable?

Did you consider the use of (potentially gated) residual skip connections (Savarese & Figueiredo 2017)? This (or a related approach) will likely help to increase/improve the expressivity of these attention layers and prevent oversmoothing by allowing for a more persistent signal from first- and second-order epistatic interactions, potentially allowing for the use of additional layers (if necessary).

Is there a significant difference in the mean level of phenotypic plasticity? It's fair to ask whether the distributions of the extent of pleiotropy is the same for convergently vs non-convergently evolved genes, but the conclusions you are reaching about a directional effect (i.e. convergent genes exhibit greater pleiotropy) are not aligned with the tests being conducted, which simply address distributional different, irrespective of direction.

How is the number of single gene translocations since divergence from their common ancestor determined? Does this assume or leverage some rooted phylogeny with an outgroup? Put another way, how can you distinguish between derived vs ancestral & conserved gene order?

Earlier in the paper you highlight that the METL-Global model tends to overfit fairly strongly to the in-distribution training data. I can't help but wonder whether the increased usage of multi-headed attention here and equal, but still relatively low dropout may be contributing to this. Did you explore the impact of increasing dropout on the ability of this model to generalize? Related, did you explore whether implementing dropout to the attention heads themselves improved the models capacity to generalize? I fully appreciate this can be tricky to nail down, but I can't help but wonder if doing so might help somewhat for out-of-distribution predictions, even if it is to the slight detriment of in-distribution prediction.

For how many variants did this occur? Was this concentrated in specific proteins in the global dataset? I'm just wondering/curious if these NaN's might be associated with some specific, biologically relevant phenomena leading to some potential acquisition bias in the simulated training data, or if these instead can be interpreted as an unbiased occurrence.

I agree with this interpretation insofar as it suggests that the sampled functional landscape for these proteins is dominated by additive effects. But, even in cases where multi-mutants have been sampled, these datasets typically only have sampled multi-mutants at a handful of sites at small hamming distances from WT.

For instance, in the case of the avGFP dataset which includes mutants up to an impressive hamming distance of 15, the sampling density of multi-mutants steeply drops off at increasing distances. Taking into account the fraction of possible multi-mutants sampled at those increasing distances, it's apparent that we've barely scratched the surface of these protein's sequence space and presumably functional landscapes.

This is all to say - these datasets (naturally, due to the nature of their generation) have very strongly biased sampling of the functional landscape, being extremely localized around a single WT protein. It may well be that these WT sequences reside on fitness peaks around which local mutational effects are largely additive.

I would be interested in seeing the regime extrapolation performance of these methods broken down by mutational distance from WT, or by the count of mutations. I can't help but suspect that the strong performance of linear models decays strongly with mutational distance. When considering that these datasets often have sampling biases towards shorter mutational distances, it seems likely that the correlation performance will be dominated/confounded by the high abundance of "easier," local multi-mutants.

How many parameters are in each model? I would suggest reporting this for each of the base models as well as their ensembles.

I would like to see figures depicting the performance (error) of these different methods for each parameter estimate, but plotted as a function of tree size (i.e. a head map of error as a function of parameter value and tree size) - I can't help but wonder if part of the pattern of increasing error rates as a function of increasing diversification parameter is simply a result of there being greater variability in simulated tree shape/size at these larger parameter values.

Additionally, I think it could be quite informative to plot a heatmap of error with speciation and extinction rates as the x and y axes - I suspect this would highlight a clear, predictable pattern, particularly with increasing error rates being characteristic of parameter combinations where both speciation and extinction rates are high, leading to high species turnover and thus greater "volatility" of diversification outcomes.

Again, I can't help but wonder/suspect that the performance of the GNN here is limited due to the fact that GraphSAGE cannot actually leverage/take advantage of edge weights/branch lengths.

Related, the performance of the GNN may be improved further if node positional information in the tree is encoded using some of the node representations implemented by pytorch geometric as node features (e.g. the Laplacian Eigenvector positional encodings - https://pytorch-geometric.readthedocs.io/en/latest/generated/torch_geometric.transforms.AddLaplacianEigenvectorPE.html#torch_geometric.transforms.AddLaplacianEigenvectorPE).

This citation should be provided for the first reference to this code base.

It's nice to see that you used AdamW over the original Adam optimizer, as it correctly implements L2 regularization. But, it might be worth considering the use of either adaptive learning rate schedulers on loss plateau to improve/generalize the ability of each of these three models to learn during the training procedure, as well as to possibly explore the impact of varying the value of weight_decay (regularization strength), as the default value for this parameter (and realistically for the learning rate) is unlikely to be the same for these three models which differ in model complexity.

Did you explore the effect/impact of ordering your model predictions on the outcomes? If so, what sort of variability in outcome did you find?

I comment on this in Appendix C, but did you explore the use of different GNN architectures? GraphSAGE does quite well for inductive learning tasks and for aggregating node information from multiple hops, but it does not take into account edge weights or edge attributes, which here could be interpreted as branch lengths - I suspect this information would be incredibly useful for a GNN model capable of leveraging this information.

Is there a particular reason why you chose to use GraphSAGE? I understand it does well in aggregating information from nodes that are multiple hops away and thus performs well in large graphs, but it is unable to take advantage of edge weights (i.e. branch lengths) which could be incredibly information rich for diversification parameter estimation. Additionally, although you capped the tree sise to 1500 tips (3000 nodes), these graphs really are not terribly large in the context of GNN applications - they are also quite sparse, as they are fully bifurcating.

Mostly just interested in your thought process here! Did you explore the use of alternative GNN convolutional layers or models?

Nit: I'd suggest writing out the units - 1.5 species / unit time. Alternatively write out mu = 1.5

I think you mean for this to be capitalized as"EvoNN" as is done elsewhere.

Is this meant to be EvoNN? Or is this a separate package/repo? If the latter, please be sure to include a link to the repo!

I think it certainly makes sense to evaluate these different model architectures for these commonly studied diversification scenarios for which existing, comparable PCMs exist. However, I do think it's going to be incredibly important to evaluate how these do under increasingly complex scenarios, such as those where diversification rates are time- and branch-heterogeneous and either depend on some other continuously (or discretely) varying trait(s), or themselves evolve according to some model of evolution (e.g. Martin et al., 2023 - https://doi.org/10.1093/sysbio/syac068).

I say this because although these diversification scenarios may require even more information-rich, larger phylogenies to accurately infer diversification parameters/make phylogenetic predictions, these are likely to be the scenarios where more classical MLE-based approaches exhibit a high degree of model inadequacy, leading these new approaches to shine.

It's great to see this recognition/appreciation for the potential of GNNs in the context of phylogenetic inference - there are a huge range of applications here!

Given the strength of the claims you are making here, I would suggest reporting these results explicitly - as is the reader has no ability to evaluate this claim.

While useful in assembling sequences for your MSA, HMM methods for MSA construction may no be, or are unlikely to be the most accurate. Given that you are interested in the impact of MSA construction on your predictions, it might be worth considering the more accurate methods typically used in phylogenetic inference.

Is there a reason why you only considered the one feature based on protein 3D structure?

Have you considered evaluating how the use of a reference-free approach (as described here) might influence your findings here? By evaluating everything with respect to the wild-type, it can be difficult to disentangle local-epistatic effects from global nonlinearity in the sequence-function relationship.

https://doi.org/10.1038/s41467-024-51895-5

While certainly of importance, I would suggest this statement be tempered slightly - I'd suggest instead emphasizing that protein stability is one of the strongest forces constraining molecular evolution.

I would suggest clarifying earlier in the results section that across all models, the feature encodings used were ESM embeddings - as is this doesn't clearly come across without first going into the methods.

I would suggest describing/defining here exactly what you mean by MSA depth - does this involve collection/use of out-of-sample natural sequences to construct the MSA? Is depth here the number of aligned sequences outside of those you are predicting fitness for?

I presume this is in reference to "Hamming distance EVmutation," plotted as the brown dotted line? This does quite well, but it's unclear from the figure with this actually is/how it incorporates the two features to make its predictions, since it's not described/indicated elsewhere in the figure.

I like this idea - but it seems like it would be quite valuable for you to explicitly assess the extent to which this is true. For instance, does the number of KDE peaks correlate with landscape ruggedness? Or the number of fitness peaks, or fitness sinks?

I might suggest adding in some of these basic exploratory tests investigating the relationships between the fitness statistics and fitness landscape measures to your work, either in the main text or supplement.

You haven't yet mentioned anything about how what ML model or what model architecture you are actually using to make fitness predictions - these details may be in the methods, but I would strongly suggest a section/paragraph briefly describing this to improve clarity!

Is there a reason why you've chosen to use "local optimum" instead of the more commonly used term of "fitness peak" in the fitness landscape literature?

Is it actually appropriate to consider these as being fundamentally distinct from each other? It seems like it would be more apt to describe them as different, but (in some cases) partially overlapping regions of the same activity landscape.

This is the first time you mention that you've calculated Ka/Ks - I would suggest first describing these analyses and the motivation behind them.

What are the red branches? Are the tree topologies for structure and sequence identical?

This is a bit awkwardly phrased - maybe reword to something like "Tandem duplications are usually highly unstable and typically disappear rapidly from their natural populations in the absence of natural selection acting to maintain them"

Given that all methods are in the supplement, I think it would be best if you specify the exact methods you're using here - are these trees produced as part of the base OrthoFinder workflow? How did you infer multiple sequence alignments?

Were these protein strain datasets pre-processed to remove any potential artificial gene duplicates. That is, identifying/removing highly redundant sequences that arise due to bioinformatic reasons, rather than actual gene duplication events? Inferences of orthogroups using tools like OrthoFinder and subsequent downstream analyses are highly impacted by the persistence of such artifacts.

Below I point out that PF+FastME seems to do well when it comes to inferring branch lengths (as indicated by low KF distances), but it struggles with topological accuracy (as indicated by RF distances), particularly for larger trees.

I wonder if you might be able recover some topological accuracy by supplementing the one-hot encoded MSA with some form of complementary, or even alternative feature representation. One potential that comes to mind is that perhaps obtaining protein feature vectors derived from a Position-Specific Scoring Matrix (PSSM) constructed from the MSA.

This would be a 20 AA x 20 AA = 400-long feature vector for each sequence (such as implemented by PSSMCOOL - https://doi.org/10.1093/biomethods/bpac008). Including this as either a substitute for the one-hot encoded MSA, or as a complementary layer may provide additionally relevant information for PF to learn, resulting in improved topological accuracy while still using FastME.

These results paired with the runtimes are really quite impressive! But the contrast between the results for KF distances as compared to RF distances are interesting, and seems like they may be worth unpacking.

In particular, it's notable that the RF distances at greater tree sizes for PF+FastME seem to converge with FastME, being greater than seen for IQTree/FastTree, with the difference increasing along with tree size.

As you say, RF is just the sum of differences in bipartitions between two trees, whereas KF considers both differences in topology and branch length. You find that PF+FastME consistently infers trees with lower or equivalent KF distances to IQTree and FastTree. But, as tree size increases, RF distances increase for PF+FastME at a high rate, exceeding those of FastTree and IQTree starting at relatively small trees (~20 tips).

Together, these results would suggest that PF+FastME estimates branch lengths well. This is maybe expected but a great thing to see, since PF is effectively trained to infer those evolutionary distances that FastME uses to infer branch lengths! However, despite accurately inferring branch-lengths, there seems to be a larger number of topological errors in the larger trees inferred by PF+FastME as compared to the other methods.

Do you have any intuition as to why this discrepancy arises? Or any thoughts on how you might modify the model/model architecture to better account for and mitigate this effect?

Seems we see the same pattern as in Fig 2 here - this time with RF distances for PF+FastME converging with FastME at slightly larger tree sizes (~50-60 tips).

We see a similar pattern as before (in fig 4) in the SelReg dataset, where KF distances show FastME and PF+FastME inferring trees with better branch lengths, but (particularly for larger trees) greater topological errors. Additionally, with increasing tree size, most methods seem to show a trajectory wherein RF distance decreases with increasing tree size, whereas PF+FastME is still increasing (for both the Cherry and SelReg datasets). I think it would be worth expanding the datasets if possible to see what happens at even larger tree sizes (e.g. 250, 500, 1000 leaves). Do these patterns continue or plateau? If so, where?

Given that you report results for both RF and KF distances above in figure 5, I'd suggest doing the same here, as it seems clear that RF and KF are capturing quite different features of topological differences, particularly for IQTree vs FastME, and PF+FastME.

I know the this is just the results section, but it would be useful if you could clarify whether you're arriving at these conclusions using statistical methods or not. As written it's a bit unclear.

Related to this, it might also be interesting to explore the use of your inferred macrosynteny as characters for phylogenetic inference! For instance, rather than using the concatenated single copy genes to infer your species tree, would use of macrosynteny lead to the same conclusions with respect to the relative prevalence of fusions/splitting events?

Is this not the same figure as in 3A, just excluding T. lapidaria an and colored according to the earthworm ancestral linkage groups?

I think in some cases these ideogram plots are somewhat redundant across figures. This is ultimately a matter of preference, but it might help to streamline things by limiting the number of figures.

Where exactly did this occur? And you point to the node? From the figure it's unclear if this occurred at the base of the clade or somewhere else.

It's a bit difficult to gain clear intuition here regarding the relative prevalence of fusion vs splitting events. Could you perhaps plot the phylogeny with nodes labeled according to the count of fusions/splits (e.g. something like 2/0, where # fusions = 2, # splits = 0)?

Where is this coming from? Are these numbers from explicit statistical methods, e.g. ancestral state reconstruction? I'd be clear about this, because as it, it reads as more of a hypothesis than an inference.

Why did you decide to use concatenation rather than a partitioned analysis, or more formal species tree models capable of leveraging multi-copy gene families such as SpeciesRax (https://doi.org/10.1093/molbev/msab365), Asteroid (https://doi.org/10.1093/bioinformatics/btac832), or Astral-Pro 2 (https://doi.org/10.1093/bioinformatics/btac620)?

I'd definitely suggest exploring the use of one of these alternatives, as the models for species trees are quite different from those for concatenation. Furthermore, in some cases concatenation can lead to highly "resolved" phylogenies with strong bootstrap support. For a more detailed discussion of the matter, I'd suggest taking a look at Liu et al 2015 (https://nyaspubs.onlinelibrary.wiley.com/doi/abs/10.1111/nyas.12747)

Using IQ-tree? I'd suggest specifying the method you chose to use here.

The order here seems a bit off, and led to one of my comments later on the that stemmed from this confusion. I would suggest ordering such that the sections describe: 1) The collection/curation and gene model annotation of the 23 genome assemblies and summaries of them. 2) the inference of orthology 3) phylogenetic inference, and 4) the remaining analyses of macrosynteny

I suggest this because currently, you describe phylogenetic analyses which use the orthogroups inferred using orthofinder, but those analyses or results are not described until after the phylogenetics. I think restructuring in this way might be a bit more intuitive for the reader.

Given the earlier comments about uncertainty in transcriptome based phylogenetics and the fact that the new phylogenetic hypotheses presented here using support those from transcriptome datasets, I might suggest either tempering the earlier statements, or describing in greater detail/specificity here exactly where there are differences in topology between the two.

As written it's a bit ambiguous what the data used for phylogenomic inference here were - I might suggest rewording slightly to something along the lines of:

"We built a maximum likelihood phylogeny of annelids using >500 single-copy genes sampled from the chromosome-level genome assemblies and newly annotated gene models"

Can you elaborate on/discuss why you believe transcriptome based phylogenomic inferences are more likely to retain uncertainty than from other genome-wide sources of data?

I ask because there will always be variable degrees of uncertainty in phylogenetic/phylogenomic inference, irrespective of the source of data used. Furthermore, (measured) uncertainty can be confounded by the methods used, source of data, and more (see Simon 2022 for a discussion: https://academic.oup.com/sysbio/article/71/4/921/5904279).

It's a bit confusing to me why you would carry out your model selection procedure without using the intended feature-set. My worry would be that certain models might perform better or worse with different feature types and that you might be missing that here. Could you elaborate on this?

Did you also consider exploring the use of regularization? I Would suggest looking into the use of L1 or L2 regularization to reduce the contribution of some of your features to mitigate the potential for overfitting.

Is there a particular reason you chose this specific PSSM profile over the alternatives?

Along these lines, it's important to make significant efforts to identify sources of bias in your training dataset and mitigate their potential impact on predictions. This is true for the training set, but it's similarly true for the test (referred to here as the validation) set - if the test set is biased or imbalanced with respect to some relevant biological feature, then the resultant prediction accuracies may not reflect true model performance.

What I would like to see explored more thoroughly here is whether there are taxonomic biases in this curated set of proteins used in training and testing of your model? If, for instance, some species/taxonomic groups are disproportionately represented in both your training and validation sets, potentially leading to elevated prediction accuracies.

This is a bit confusing - common terminology to use here would be to refer to this as the training dataset (subsequently subdivided into the K-folds for cross-validation), and the latter as the test set, rather than the validation set as you've done here.

Can you elaborate on this some? As is it's not clear whether all of the listed studies above have evaluated their methods on unreviewed data, and for those that have, what the sources of the unreviewed data are.

It surprises me that you've chosen to use parsimony to infer ancestral character states when the rest of this paper and the motivation of the CAnDI software itself is clearly motivated by explicit, formal models of evolution.

If you were to infer ancestral character states using a Markov model of discrete states, are there major differences in ancestral state reconstructions as compared to the parsimony reconstructions?

The trait distribution is quite simple (not a lot of within clade trait heterogeneity), but at the very least I might suggest conducting formal model comparison among (for instance) an equal-rates, symmetric / asymmetric rates, and all-rates-different Mk models and reporting the inferred character states under these model and compare them to your parsimony reconstructions. If they are concordant, that's fantastic! But if not, it would be worth commenting on how this might impact your interpretations.

Can you provide some more specific explanation as to what this script is doing under the hood?

For instance, what characteristics lead you to classify gene copies as lower quality duplicates? How could these decisions impact or influence your inferences surrounding conflict, and could removal of legitimate but divergent duplicate gene copies bias your results in some way?

I'm not entirely clear what you mean to say with this - could you perhaps use other language to help clarify?

Do you mean to communicate that the method allows users to identify specific instances or concerted patterns of conflict that are coincident with or predictive of trait evolution?

Given that the dates inferred from your time calibration are of such critical importance to the conclusions made herein, I might strongly suggest using more robust methods that can take advantage of your fantastic UCE dataset. In particular, I would recommend looking into either using BEAST or MCMCTree.

Which dates are you reporting here and throughout? Presumably these are dates from only one of the two methods (Chronos and LSD2).

What were the times inferred from each method respectively?

I understand that the detailed methods are in the supplement, but you really should at least state which methods you are using, otherwise this statement is not particularly informative.

One thing I did notice when reading the supplement - it doesn't really make sense to use the IQ-TREE implementation of LSD2 here if the only constraints you have are for the root age. This method is primarily intended for use in phylodynamic contexts wherein you have tip-dates for most or all of your rapidly evolving pathogen samples. At its core, the method dates the tree by regressing the tip-to-root distance against sampling date - if you only have a date for the root, the method is unlikely to be accurate.

From the LSD2 github page (https://github.com/tothuhien/lsd2): "The input date file should contain the date of most of the tips and possiblly some internal nodes if known."

This plot is a bit misleading, as diversification events plotted here include both geographic subdivision within species (e.g. Entacmaea) as well as diversification among species. Furthermore, without seeing the clownfish phylogeny, it's quite difficult to correspond the tempo of their diversification to the tempo of diversification in anenomes as seen here.

I would strongly suggest plotting the clownfish phylogeny alongside this one as "co-phylogenies", such as implemented within phytools: http://blog.phytools.org/2016/08/empirical-co-phylogenetic-cophylo.html

If you have information about which specific host species each clownfish is paired with, this would be a useful way of more intuitively visualizing whether the diversification of mutualistic species pairs truly is coincident.

This statement really is quite strong, however I'm not convinced that it's warranted simply given a coincident interval in which extant species of both groups diversified.

For instance, given the close ecological interactions among the anemones and clownfishes, they certainly share many ecological conditions. From this, it seems clear that although one hypothesis to explain the pattern of co-diversification is that of coevolution, an alternative hypothesis is that each diversification in each group was driven by some shared feature of their life history, independent of the interactions between anemones and clownfishes.

I might suggest including a description of alternative hypotheses, and the extent to which your data and analyses support/refute them.

Again, I'd strongly recommend reanalyzing these data (particularly 3B) using phylogenetic comparative methods! Although some taxonomic groups seem to exhibit a strong trend (e.g. flowing plants) others do not.

This is where I would strongly suggest reanalyzing these data using phylogenetic comparative methods - you may be surprised to find that previously non-significant relationships become significant after using the more appropriate methods!

Are you referring to genes here? or the genome? In other words, do you mean to say that less than 50% of genes in most multicellular life forms are coding? Or do you perhaps mean to say that less than 50% of the genome is coding? If the latter (which is what I suspect you mean to say), I would revise throughout for clarity/accuracy.

This figure demonstrates clearly the strong effect of shared evolutonary history (i.e phylogeny) on the relationships recovered herein. Another potentially interesting analysis you could consider here is a phylogenetic analysis of covariance (phylANCOVA - originally published by Fuentes-G et al., 2016 - https://doi.org/10.1086/688917).

This test is effectively equivalent to a normal ANCOVA in that it formally tests whether different groups (here, taxonomic groups) are best explained by different regression models. In other words whether groups share slopes, intercepts, or both. The difference is that the test can account for evolutionary non-independence using a phylogeny - this test could be conducted using the a phylogeny comprised of the full suite of species for which you have these data.

This could be a nice way to formally assess which groups should actually be modeled separately, and assessing whether differences among groups are statistically significant. Currently, difference among groups is implicitly assumed, by modeling each using their own model.

Minor typo suggestion - throughout, "taxonomical" should be rewritten as "taxonomic".

To conduct your regressions in a manner that is statistically robust to evolutionary non-independence within each taxonomic group, I suggest you conduct phylogenetic least squares (PGLS). This can be conducted in a very straightforward manner using R (for instance, as described here: http://www.phytools.org/Cordoba2017/ex/4/PGLS.html). All you will need is: 1. The per-species measurements (i.e. genome size, gene content, length of coding sequence, etc), and 2. a time-calibrated phylogeny (branch lengths are in units of time) that contains the species for which you have these data.

Since you have measurements for so many species (a fantastic dataset!), I might suggest simply curating your list of species and obtaining a phylogeny from timetree.org. These phylogenies will not necessarily be the most accurate, as they are effectively aggregated from the published literature, but they should be sufficient for your purposes. Although these phylogenies my not include all species in your dataset, they may include closely related, representative species that you can use in their stead.

Even if this approach requires the removal of some species in your data, use of phylogenetic comparative methods like PGLS are critical for the statistically robust analysis of these data, and thus interpretations reached herein. Without conducting these improved analyses, it unfortunately will be impossible to evaluate the accuracy of the results presented here, particularly with respect to parameter estimates, and differences in the role of alternative splicing in driving the evolution of genome complexity.

I really do like the idea of the metric you've come up with here, but my major concern is that in these correlative tests you have not accounted for the confounding effect of evolutionary non-independence between species (Felsenstein 1985: https://www.jstor.org/stable/2461605).

That is, because species share evolutionary histories, more closely related species will be more similar in their phenotypes and relevant to this study, in their genome content. This induces, in effect, statistical non-independence of species - phylogenetic pseudoreplication - that must be accounted for formally in your statistical analyses.

Failure to do so can lead to (for instance), biased or incorrect parameter estimates (both intercept and slope) within models (i.e. within taxonomic group), and can confound or mislead interpretations among groups (i.e. among models fit for each taxonomic group). Fortunately, the solution is relatively straightforward - I've detailed some fixes in the sections below.

I would suggest avoiding language such as this, since it alludes the the inaccurate "Scala Naturae"-type ladderized view of evolution, progressing from "low" complexity organisms towards "higher" complexity multicellular organisms.

Maybe specify "of extant genes" to clarify how these two processes are not inclusive of the formation of ALL new transcripts, including those for genes not yet contained within the genome of a given lineage. That is, horizontal gene transfer and de novo gene birth are also two common mechanisms by which new genes (and thus transcript diversity) arises.

It's very relevant whether a species tree vs "specimen tree" was used here, as the inclusion of specimen here could change interpretation of rates quite significantly

Okay, yes, so I'd check to see whether the parameter estimates under a model without the jump parameter are consistent with these results, or if the J model has led to very different biological conclusions.

I'll reiterate again (and hopefully I'm mistaken!) that if this tree includes both distinct species and multiple specimens from the same species, application of ancestral range (or state) reconstruction is not correct. These methods should only be applied to a species or a uniform taxonomic rank/biological unit

Am I correct in understanding that in this hybrid-ASTRAL tree, you have collapsed specimen into distinct species? If so, how many species in total?

If not, this poses problems for the ancestral state and ancestral range reconstruction, as each should only be applied to a single taxonomic rank (i.e. species), rather than multiple (i.e. species and individuals sampled from different populations)

This is not reported in the methods section, but should be.

Additionally, if I'm understanding correctly, is each tip in this tree a specimen, not a distinct species? If so, it does not make sense to use a lineage through time plot here unless you have retained only a single representative per species.

Related to my other comment, if the best-fit model was one that included the "jump" parameter "+J", this would be consistent with a demonstrated pathological favoring of such models, independent of whether such a model is biologically sensible.

See Ree and Sanmartín 2018 for more: https://doi.org/10.1111/jbi.13173

Which models did you fit? These details - both the models you fitted, and their relative fits to the data - really must be reported here in the methods. Unless I'm missing something, they don't seem to be reported in the supplement either.

I'd imagine these should have bayesian credible intervals as will, since they were inferred using MrBayes, no? It would be good to report these throughout.

"nearly an"

Also, was this inferred to have been a single origin?

Gene loss and gene gains are just the mechanisms by which (gene) copy number variation arises "Maybe reword this passage as: "but can also arise through distinct molecular paths and by distinct evolutionary mechanisms, such as copy number variation that arises via either gene losses or gene gains (e.g., horizontal gene transfer or HGT)."

I think it would be worth citing this recent excellent paper that synthesizes the many different perspectives on convergence vs parallel evolution, etc. The framework they outline for distinguishing these different processes would be useful to borrow from throughout, particularly as you begin to delve into the hierarchical nature of convergence (described nicely in Rosenblum et al., 2014).

I'd suggest just clarifying up front that convergence specifically refers to the case where similar traits independently evolve from distinct ancestral states, as is commonly the case in distantly related species. This differs from parallel evolution, in which similar traits evolve, but from the same ancestral state - this can still occur among distantly related species, but more commonly manifests at the population level, or among more closely related species.

This could be as simple as: "Convergent evolution, the repeated evolution of similar traits from distinct ancestral states, is ubiquitous in nature, commonly occurs among distantly related taxa, and has been..."

I'm realizing as I read through the methods that it is not described here how πS was actually calculated. This is so crucially important for the paper, so these details really must be included in detail. For instance, reading the methods it's unclear to me whether πS was calculated in sliding windows across the genome, using the core genes, or something else. The details may be outlined elsewhere in the manuscript, but the really must be here in the methods. It's quite difficult to evaluate the results otherwise.

I would suggest being more specific about what you mean here - which specific parameters are you referring to as "adaptation rates"?

Why were these methods applied to the concatenated genes? Would it not be more informative to fit these models on a gene-by-gene basis so as to gain more insight into how the DFE varies within and among species? Alternatively, if the issue is a matter of statistical power, could this be instead be done on the concatenations of the three datasets described above? (non-secreted, secreted, and effectors)

I would be curious to see whether Ne or the recombination rate vary in some predictable way among the species corresponding to each model.

It seems to me that accurate estimates of Ne will be critical to your interpretations throughout the manuscript, will it not? If so, and given the incredible dataset you're working with here, I can't help but feel that using πs here rather than actual estimates of Ne is surprising.

At the very least, I might suggest using a method such as GONE (https://academic.oup.com/mbe/article-abstract/37/12/3642/5869049) to estimate Ne for each species - at the very least to demonstrate that the two approaches leads to comparable or equivalent solutions. The method performs quite well for estimating contemporary effective population sizes from modest numbers of samples. a comparison of this and other related methods (including those that better infer historical pop sizes) can be found here: https://academic.oup.com/mbe/article-abstract/37/12/3642/5869049

Horizonal transfer (including introgression following hybridization) is another potential explanation - this would be a good place to discuss this point.

This needs to be explained further. How exactly was this done? Did you convert the inferred species tree to be ultrametric and then change the root depth to 87Myr? If so, this is not an appropriate method for time-calibration.

At a minimum, I would suggest either using this root age and r8s as implemented in geiger (https://rdrr.io/cran/geiger/man/r8s.phylo.html), or PATHd8 (https://www2.math.su.se/PATHd8/). Alternatively you could use another species tree time calibrated with other more intensive methods that includes some proportion of these species, and use the congruificaton method (again implemented in geiger (https://rdrr.io/cran/geiger/man/congruify.phylo.html) which uses PATHd8 or treePL under the hood.

If you've used one of these methods to time calibrate your tree, be sure to describe that here.

How sensitive might your results be to the inferred species tree topology? Numerous efficient methods exist that are capable of inferring species trees from multiple copy gene families, including asteroid (https://github.com/BenoitMorel/Asteroid) and ASTRAL-Pro2 (https://github.com/chaoszhang/ASTER/blob/master/tutorial/astral-pro.md).

Additionally, SpeciesRax (https://github.com/BenoitMorel/GeneRax/wiki/SpeciesRax) is capable of inferring a rooted species tree using multi-copy gene family trees (in this respect similar to the STAG algorithm implemented in OrthoFinder) under a model of gene duplication, transfer, and loss.

I would suggest applying one of these methods to assess how robust your species tree topology is to choice of methods, as the inferred gene family evolutionary dynamics will naturally be sensitive to this.

Could you briefly explain why this approach was chosen rather than species/gene tree reconciliation methods like ALE (https://github.com/ssolo/ALE) or GeneRax (https://github.com/BenoitMorel/GeneRax)?

Because CAFE only models gene duplications and losses, it's unable to account for or infer potential horizontal gene transfer among your focal species, a process that could (and likely does) explain some proportion of the variation in gene copy number observed in the analyzed gene families.

I would suggest briefly addressing these points, first in the methods section, and subsequently in the discussion. HGT is a plausible alternative hypothesis to gene duplication, and one that is sensible in the context of gene family-diversity associated increases in a species' capacity for hostplant range expansion.

This makes good sense, though I might suggest elaborating on what sorts of biases high variances might introduce (e.g. 'cascading' influence on inferred expansions on branches neighboring unsusually high copy numbers?).

Additionally, how exactly did you filter - did you calculate variance among species (or std. deviation for instance) and then filter based on some empirical threshold?

How exceptional is this count of shared gene families among the two symbiotic lineages? Does it exceed what is shared between only (for instance) Ev and Po, or S1/2 and Po?

I feel that convergence of the gene families shared between the two symbiotic clades is an alternative hypthothesis that should be tested, rather than assumed. An alternative explanation is that these are genes that arose once within S1, and were subsequently lost within Ev.

Again, I would just be cautious in the assumption of convergent evolution here, as this is an alternative hypotheses that has not yet been explicitly tested.

This sub-plot would again support an interpretation in which (at least some/many) genes related to a symbiotic life history were gained once in S1, and subsequently lost in Ev.

Okay, this addressed my comment above, but I would say that this supports an interpretation that the large numbers of genes shared between the symbiotic clades arose once in S1, but were subsequently lost in Ev, being retained in S2.

I think this is a case in which using explicit phylogenetic comparative methods might be useful. If all species were included in the species tree (rather than collapsed by clade as discussed above), and intron lengths were mapped onto the terminal branches of that phylogeny, ancestral character reconstruction could be performed to gain more explicit intuition as to whether independent expansions or contractions occurred, and how many times.

This could feasibly be done both using continuous measures of intron length, or alternatively by binning introns into "small" and "large" size categories (if intron length follows a bi-modal distribution for example), and using a discrete model of character transition to explicitly obtain counts of the number of transitions.

Is there a reason why S1/S2 were collapsed, rather than showing all sampled Suessiales species the species tree? It seems this is a loss of information that would be nice to see in this early figure.

In particular, since Symbiodinium is not exclusively symbiotic, it would be informative to see how this life history is distributed within the group phylogenetically. Different phylogenetic histories and distributions of symbiosis could lead to alternative interpretations regarding the evolution of the life history, and consequently of the evolution of parasitism-associated genome size reduction.

Given that amino acid sequences were back translated to nucleotide alignments, it would be worth assessing to what extent the tests of relaxed selection may be sensitive to/impacted by alternative translated alignments. That is, each amino acid sequence has multiple alternative nucleotide alignments.

I'm not certain how pal2nal works exactly, but from what I can tell it appears to generate a single translated alignment. Depending on how it "chooses" which alignment to generate, this could artificially inflate (or deflate) signal of relaxed or positive selection, as it may erode the very signals these tests use to fit model parameters.

Again, just wanted to give a heads up that this seems like this repository may still be private!

This is not terribly surprising to me, honestly, particularly given that it's known a priori that the two lineages are closely related and that hybridization is not uncommon. It may also be worth mentioning/emphasizing here that the phylogeny was inferred using concatenated SNPs. Whether using full sequences or SNP datasets, concatenation can often lead to inflated topological support.

Have you considered using an approach like SVDquartets to infer a complementary phylogeny to the one inferred with RAxML? The approach is consistent with the multi-species coalescent, and treats sites as independent, inferring quartets at each position before amalgamating them into a full tree. Individuals could be pooled by sampling locality/species, enabling you to still infer a population-level tree, doing so under an independent phylogenetic model.

I'm super intrigued by the use/utility of simulation here. I know that these simulations were conducted using msprime, but I can't help but wonder about the inclusion of varying intensities of positive selection on loci in these simulations. This of course would increase the complexity of parameter space significantly with regards to potential simulations, but seems explicitly tied to the hypotheses being tested here.

More generally, it seems to me that this framework could lend itself nicely to the application of machine or deep learning approaches for assignment of genomic windows to alternative evolutionary histories (e.g. migration, selection), as well as potentially even parameter inference, such as through the use of a CNN. Obviously this is outside of the scope of the current paper, but I'm just curious as to whether this is a thought/space you all have explored?

First off, I can't not complement the "Twisst & Tern" pun... Excellent stuff.

But more specifically, I just want to say I think the use of the ternary plot in combination with both empirical and simulated data to explore and quantify genealogical asymmetries is so clever. I doubt I'll be alone in saying that I think this type of approach has immense potential, not only for hypothesis testing, but even for parameter testing.

Interesting - this seems like a perfect complement for the analysis of Cocoonases later on! Would be very interesting to apply RELAX/aBSREL in the hypthesis-testing framework (i.e. foreground/background - discussed in later comment) to this gene family to test the hypothesis that this gene-family contraction in H. aoede is paired with an increased intensity of positive selection as compared to other Heliconius butterflies.

What sre the members of the Hemocyanin superfamily that contracted within H. aoede? Are they the same (or do they include) the two Hemocyanins that were found to have previously undergone expansions?

Under what model specifically was ancestral state reconstruction (ASR) conducted? The ML method in phytools can conduct ASR under three models (Brownian Motion, Ornstein–Uhlenbeck, and Early Burst). If not already done so, I might suggest fitting these alternative models of trait evolution prior to conducting ASR under the best-fit model.

Given the taxonomic scale, it might also be worth considering including fitting some multi-regime BM/OU models as implemented in OUwie: https://doi.org/10.1111/j.1558-5646.2012.01619.x.

If I'm reading this correctly, it seems that for each analysis, all branches were tested for signatures of selection (with the exception of SLAC which assumes selection is constant across the entire tree)? Although RELAX and aBSREL can be used to test each branch, each are significantly more powerful when used in a hypothesis-testing framework, assigning a subset of branches to the "background" (control) and the remainder to the "foreground" (treatment), and estimating parameters within each test-set. The per-branch tests are generally considered to be most useful in exploratory use-cases.

I think this use-case scenario is ideally suited to the hypothesis testing framework of these two models. Specifically, to test the hypothesis that selection on cocoonases intensified (or was relaxed) in H. aoede as compared to the rest of the Heliconius butterflies, you could treat H. aoede as the foreground, and the rest of Heliconius as the background for both RELAX and aBSREL. It's unclear whether cocoonases in Eueides or other species within Heliconiini should be experiencing the same/similar selection regime as in H. aoede, so I would suggest excluding those species from this specific analysis.

How might this approach to identify orthologs, which is tailored to the identification of single-copy orthologs, impact downstream inferences about patterns of gene-family evolution? Put another way, if the approach described here is more prone to recovering scOGs (or low-copy number OGs), then it seems likely that fewer multi-copy number (particularly large/high-copy number) OGs will be identified as a result. If so, the rate estimates of gene-family expansion/contraction using CAFE seem likely to be downwardly/upwardly biased respectively, as the majority of gene-families under consideration are particularly conserved with respect to copy-number.

I think this possibility would be worth exploring/addressing. For instance, if the implementation of broccoli is relaxed so as to be less tailored to the recovery of scOGs, does the shape of the gene-family expansion/contraction distribution fundamentally change?

Really exciting work, and an immense amount of work - congratulations! This was a fun read!

It seems that towards the conclusion of the paper, you outline how your data does not appear to be consistent with one of the prevailing hypothesized mechanisms for the digestion of pollen in Heliconius. I can't help but wonder if, using the data you present here, you might be able to point towards individual gene-families/orthogroups that are consistent with the evolution of this key-innovation?

What does the distribution of mean per-species gene-count look like across all orthologs? Is there a small handful of very large orthogroups that are observed in all species and comprise the majority of the gene-counts shown here?

What was the estimated ω for H. aoede at these gene-families as compared to those distributions estimated for Eueides and Heliconius? If involved in the evolution of pollen-feeding, I think that we would expect ω to be more similar between H. aoede and Eueides than between H. aoede and Heliconius?

These two relationships should probably be assessed in a phylogenetic context, accounting for shared evolutionary history and non-independence among species (e.g. using phylogenetic least-squares regression - pgls). There is clear phylogenetic structuring in these two plots. For instance, Dione/Agraulis appear to have lower TE content for their genome sizes (shallower slope and intercept), whereas Erato seems to accumulate TEs more rapidly for their corresponding genome size as compared to other clades. Raw points and the original regression slope certainly provide taxonomic intuition, but to summarize the relationship across the entire clade, a regression slope as inferred by PGLS (for instance) should be reported.

It's difficult to contextualize these results without access to the supplementary materials (at the time I'm reading this), but it's unclear to me why the choice was made to subset introns into those that are short and long (as in Fig. 4A). I'm particularly wondering about the choice to divide introns as short/long at the median, which doesn't appear to sort them according to a natural break in the length-distribution. Dividing introns in this way, and then summing for each species and each 'TE size-class' runs the possible risk of exaggerating the observed difference in intron length and TE content.

Additionally, the relationship between TE content and total intron length (regardless of how these traits are defined) should be analyzed in a phylogenetic context (e.g. using phylogenetic least-squares regression - pgls).

Could removal of in-paralogs influence/impact the signatures of selection that you recover using aBSREL and RELAX? It might be worth assessing whether these results are robust to the inclusion of in-paralogs.

For instance, neofunctionalization of in-paralogs following duplication is likely to manifest signatures of positive selection (Fernández et al., 2020: https://doi.org/10.1093/molbev/msaa110, Wang et al., 2015: https://doi.org/10.1007/s11103-015-0285-2). Alternatively, pseudogenization of in-paralogs is expected to lead to a relaxation of selection on the duplicated gene-copy (e.g. Emerling 2018: https://doi.org/10.1016/j.ympev.2017.09.016; Calderoni et al., 2016: https://doi.org/10.1038/hdy.2016.59).

I'm curious - If you were to simulate sequences/trees under a model of mass-extinction and early burst but constant substitution rates, how often (if at all!) would Janus infer substitution rate shifts near the extinction event?

Could the extinction-driven bottleneck and heightened species turnover in the subsequent radiation amplify apparent differences in nucleotide composition even with molecular rates remaining constant? If so, could it be that alternative diversification dynamics/models impact the inference of substitution rate shifts?

Is there any chance that this pattern could, at least in part be attributed to the relative paucity of branches in the tree prior to the K-Pg boundary?

In other words, could the identification of molecular rate shifts by Janus be impacted by "sample size" differences of branches prior to and following mass extinctions as in this example?

Super cool - this is partly where/why the use of phylogenetic path analysis may be particularly informative - could these shifts in evolutionary optima have been associated with a shift in causal relationships among biological/life history features? Maybe it's possible to infer, at least in part, the causal mechanisms behind these shifts in metabolic allometry?

I'm not necessarily suggesting you do this, since this paper already represents an impressive effort, but I wonder whether the use of phylogenetic path analysis (https://doi.org/10.7717%2Fpeerj.4718) might be informative here?

Within each inferred molecular model shift, you could test alternative causal models depicting the relationships between rates of molecular evolution and these natural history traits, as you might anticipate that many of these covary, themselves being causally related to one another. Inference of shifts in causal relationships among natural history traits in each clade that experienced shifts in rates of molecular evolution could then provide some intuition as to why!

I'm probably sounding like a broken record now (a good thing - this is super interesting to me, and it's exciting to be nearing a place where we can address this!), but I think analysis of alternative causal models here with phylogenetic path analysis would be an awesome way towards not only generating this intuition, but also explicitly testing these hypotheses!

How feasible would it be to infer these deviations from equilibrium base frequencies from simulated data? For instance, simulating trees/sequences under constant substitution rates (e.g. using AliSim - https://doi.org/10.1093/molbev/msac092) but under a model of mass (selective) extinction and subsequent early burst of diversification (e.g. TreeSim - https://doi.org/10.1093/sysbio/syr029)?

I ask because I wonder whether simulating under the hypothesized diversification model and constant rates of molecular evolution could provide an estimate of the false-positive rate with which Janus might infer these shifts and whether false positives tend to be concentrated around the time of the simulated mass extinction?

I realize this is probably a big ask though - the paper is already super impressive! I think this analysis could just be particularly compelling evidence if in fact this concentration of events is greater than expected.

Would be exciting to explicitly test this!

That's one possibility - though it might be worth addressing/discussing the alternative hypothesis - that the pattern of mass (and presumably non-random with respect to body size) extinction and subsequent radiation (i.e. the early burst) could lead to surviving lineages exhibiting greater disparity in rates of molecular evolution and other life history traits.

As I alluded to in my comment on figure 1, it might be worth (if feasible) explicitly assessing the probability of observing such a concentrated distribution of shifts near the time of a mass extinction event by simulating under a constant rate model of nucleotide substitution, and a mass extinction/early burst of speciation/turnover.

I have a couple thoughts about this - is there a reason why a combination of fastme and RAXML was used as opposed to full likelihood inference using RAXML or IQ-TREE?

I ask because both RAxML and IQ-tree include methods (that i would suggest using here) for acquisition bias correction when using SNPs that could improve inference of both topology and branch-lengths.

Another entirely different strategy to infer a species tree in a computationally efficient way using SNP data under the multispecies coalescent would be to use SVquartets (Chifman and Kubatko 2014). This potentially well suited for your dataset as it can leverage a large number of SNPs effectively, and even robustly handles linked sites. Caveat is that it cannot infer branch-lengths, but this might not be a concern for your use case.

I may have missed it, but I don't believe it is mentioned in the methods what software was used for phylogenetic inference of the recovered UCEs? If my intuition is correct, I would guess IQ-Tree, the same as for mitochondrial gene? It might be worth clarifying here.

Related, was a single concatenated phylogeny used? Or did you infer gene trees as well? IQ-Tree is capable of easily inferring both (gene trees and concatenated species tree - see here: http://www.iqtree.org/doc/Concordance-Factor). This would also enable to you infer multiple concordance factors which could help to summarize genealogical concordance, and potentially even quantify per-site discordance using your concatenated SNP dataset?

Could you clarify why 15?

I'm confused though how coexistence without hybridization of closely related species of the same ecomorph suggests character displacement?

Maybe this is clarified later (e.g. some other feature of morphology/ecology other than coloration differs within ecomorph), but if so, it should be explained/mentioned earlier in the paragraph.

For instance, do spiders vary morphologically within and among ecomorphs equally? greater? less than?

Given that coloration is the only phentype explored here it might be worth providing some additional context to this point - as written it suggests that the paper might later investigate character displacement, though analyses are primarily genetic.

Could you elaborate on why you chose to use this tree versus the UCE likelihood tree? Are results from Dsuite robust to the choice of species tree used?

What is the extent/degree of within versus among ecomorph morphological variation? Is ecology and coloration the only characteristic that reliably distinguishes them?

This is a very specific number to use as a filter. Was this determined based on patterns of LD-decay inferred from the steps above? I might provide some additional explanation as to why this was used as the LD-filter.

If i'm understanding the description of the UCE tree correctly, would it be correct to say that there is one additional mitonuclear discordance when considering the UCE tree? (The placement of T. anuenue)

Oh, I understand now! It might be worth adding this explanation in the methods? I would also suggest including a plot of LogLik for each value of K - it's difficult to evaluate how strong of a difference there is in model supports without it.

However, it is somewhat surprising that that these data would suggest there are only 2-3 clusters, when there are 15 described species - this seems like quite the discrepancy.

Which nuclear tree? I see below that there is some discordance in topologies inferred using the two nuclear datasets (SNPs with fastme, UCEs with IQ-Tree), so it would be worth specifying here.

Was this quantified somehow? I might have missed it, but I don't see this described in the methods.

Just a suggestion, but it might simplify things for you if you had a section prior to any use of ANGSD describing the commonly used parameter settings across these analyses?

It might be really interesting to see the distribution of D-statistics pooled into both within vs among species and ecomorph comparisons.

I say this because it could be really informative to compare/rank the evidence for hybridization across levels - within/among species, as compared to within/among ecomorphs.

Could you clarify what you mean by this? It's difficult to intuit what this discordance looks like without the supplement being available. Regardless, it's nice to see that the discordance between the two topologies seems modest -

Could this not be said about gene content evolution as well? Patterns and dynamics of gene family gain/loss is certainly expected to be a heterogeneous process - does the reversible binary substitution model allow for site heterogeny?

This is the same type of gene-content encoding used by Ryan et al. 2013 as their second dataset, which coincidentally supported the hypothesis of Ctenophores as sister.

I'm curious as to why an alternative recoding of these gene families was not considered - gene copy number. From the analyses conducted here to infer orthogroups using OrthoFinder, this encoding would be fairly straightforward to accomplish, and could perhaps be even more information rich than presence absence. Additionally this would be an entirely distinct way of encoding these data as compared to the approach employed by Ryan et al., (2013).

Additionally, there exist now several methods that are capable of inferring species trees from multi-copy gene family trees, which are by design more robust to, and even improved by the existence of paralogs (e.g. Asteroid -https://doi.org/10.1093/bioinformatics/btac832, SpeciesRax - https://doi.org/10.1093/molbev/msab365, ASTRAL-Pro 2 - https://doi.org/10.1093/bioinformatics/btac620).

This is an impressive paper containing a significant amount of work, and I love the amount of effort focused on exploring how sensitive the results are to methodological approaches and parameter specifications. In particular, I'm delighted to see the exploration of how the choice of inflation parameter impacts downstream phylogenetic inference - given that this parameter has profound influence on all downstream analyses, I think all studies using these types of approaches should be similarly cautious/considerate.

In general, I feel that the choice to encode 'gene content' as gene presence makes distinguishing these results from past efforts, particularly that of Pett at al., 2019, more challenging. That is, gene content is encoded in the same way across the two studies - and both the present study and that of Pett et al., come to the same conclusion that the Porifera-sister hypothesis is most strongly supported. Together, it makes these new phylogenetic analyses of gene content less compelling as corroborating evidence.

I feel that an alternative encoding of gene content - as gene copy number as an ordered-discrete character e.g. (0,1,2,3, etc) - would have been valuable to explore. This type of encoding would be more information rich, and would be inherently distinct from past efforts. Alternatively, the authors could have used a growing number of methods that are capable of inferring species trees from multi-copy gene families (see later comment). I would hope to see some discussion of these alternatives, though i suspect implementation of either would be non-trivial.

I am delighted to see such a thorough exploration of the sensitivity of phylogenomic inference to these parameters, as defaults are typically used without consideration but can, as you say, have profound influence on downstream outcomes.

Unclear a priori why this might not always be expected to be the case, even when encoding data in other ways. Fundamentally, these are still evolutionary models.

I think a graphical depiction of this may be useful. As written, it's unclear how you are encoding orthogroups within homogroups.

Based on your description below, it appears to me that the only difference is in the implementation of MCL clustering to sequence similarity, with orthofinder2 used to infer orthogroups (i.e. normalizing similarity scores to account for sequence length), and homomcl for homogroups). Your statement that orthogroups are nested within protein families doesn't appear to be consistent with this? Am I mistaken in my interpretation?

Given that you are using transcriptome assemblies and presumaby amino acid sequences since you're using the "protein" method, it is important that some degree of protein redundancy reduction takes place prior to analysis with BUSCO (i.e. using CDHIT or retaining only the longest protein per assembled gene). Use of transdecoders '--single_best_orf' will still retain alternative isoforms which will lead the count of "duplicates" by busco to be inflated.

The supplementary material hosted on biorxiv doesn not seem to include these scripts/parameters.

I'm also curious: on what basis were were parameter combinations judged to be more optimal than another? I like that you explored how these parameters impacted phylogenetic inference, but it's unclear what the actual optimality criteria were.

Based on the supplement it seems you defined these models as a set of three ordered, discrete states (i.e. A <-> B <-> C). Did you consider/did you fit another set of models where there could be a two-state jump (i.e. A <-> C)? These might be worth fitting/considering, as it's not inconceivable that such transitions might occur - likewise, this additional type of transition may be a bit less sensitive under scenarios of incomplete sampling with respect to the focal trait (i.e. in cases where transitions did in fact occur from A -> B -> C, but the intervening species in character state B were not sampled (or trait data is missing).

Are these phylogenetic trees the species trees? or all gene family trees?

If the species tree, were these single-copy ortholog multiple sequence alignments for each orthogroup concatenated for inference with IQ-TREE, or using multiple partitions? Based on parameters, this would suggest the former (concatenation), but I would suggest being explicit in the methods.

So transitions between A & C were included in these Mk models? I suggest including these specifics in the methods.

Here again, it will have been quite important that protein redundancy (i.e. retaining only one protein per assembled gene) has been reduced, as even these methods to prune gene families down to single copy orthologs will be more prone to inclusion of paralogs when alternative isoforms persist in a transcriptome.

When running orthofinder (or any clustering-based ortholog identification method) on transcriptome datasets, performance will be contingent upon the filtering strategy used.

I would suggest reporting all standard BUSCO estimates: Complete single copy, Complete duplicated, fragmented, and missing. Partitioning out into this resolution provides a more complete understanding of transcriptome/proteome completeness/quality.

As I describe in the methods, it is important to reduce protein redundancy (using CDHIT, retaining the longest protein per assembled gene) prior to running busco - if you don't the number of duplicated buscos will be artificially inflated.

Was a second round of orthogroup inference conducted by including an additional set of mammalian transcriptomes?

The methods used to perform this analysis are not described in the methods section - these should be included/described. Without inclusion of these methods, it's impossible to meaningfully understand/interpret these results.

I strongly suggest providing specific details of what is involved in this pipeline. What method is used for assembly? Are any steps conducted for reducing the redundancy in these assemblies (e.g. using CD hit to remove transcripts that are highly similar in sequence)? As is, the reader needs to go to the article you cite here, and then follow a hyperlink: http://planmine.mpi-cbg.de/planmine/PlanMine_Help.html#assembly

This is purported to include a description of the assembly pipeline, but is in fact a dead hyperlink.

As is, there is no way to interpret the transcriptome assemblies produced here.

I strongly suggest providing specific details of what is involved in this pipeline. What method is used for assembly? Are any steps conducted for reducing the redundancy in these assemblies (e.g. using CD hit to remove transcripts that are highly similar in sequence)? As is, the reader needs to go to the article you cite here, and then follow a hyperlink: http://planmine.mpi-cbg.de/planmine/PlanMine_Help.html#assembly

This is purported to include a description of the assembly pipeline, but is in fact a dead hyperlink.

As is, there is no way to interpret the transcriptome assemblies produced here.

Was a second round of orthogroup inference conducted by including an additional set of mammalian transcriptomes?

The methods used to perform this analysis are not described in the methods section - these should be included/described. Without inclusion of these methods, it's impossible to meaningfully understand/interpret these results.

So transitions between A & C were included in these Mk models? I suggest including these specifics in the methods.

Here again, it will have been quite important that protein redundancy (i.e. retaining only one protein per assembled gene) has been reduced, as even these methods to prune gene families down to single copy orthologs will be more prone to inclusion of paralogs when alternative isoforms persist in a transcriptome.

When running orthofinder (or any clustering-based ortholog identification method) on transcriptome datasets, performance will be contingent upon the filtering strategy used.

I would suggest reporting all standard BUSCO estimates: Complete single copy, Complete duplicated, fragmented, and missing. Partitioning out into this resolution provides a more complete understanding of transcriptome/proteome completeness/quality.

As I describe in the methods, it is important to reduce protein redundancy (using CDHIT, retaining the longest protein per assembled gene) prior to running busco - if you don't the number of duplicated buscos will be artificially inflated.

Based on the supplement it seems you defined these models as a set of three ordered, discrete states (i.e. A <-> B <-> C). Did you consider/did you fit another set of models where there could be a two-state jump (i.e. A <-> C)? These might be worth fitting/considering, as it's not inconceivable that such transitions might occur - likewise, this additional type of transition may be a bit less sensitive under scenarios of incomplete sampling with respect to the focal trait (i.e. in cases where transitions did in fact occur from A -> B -> C, but the intervening species in character state B were not sampled (or trait data is missing).

Are these phylogenetic trees the species trees? or all gene family trees?

If the species tree, were these single-copy ortholog multiple sequence alignments for each orthogroup concatenated for inference with IQ-TREE, or using multiple partitions? Based on parameters, this would suggest the former (concatenation), but I would suggest being explicit in the methods.

The supplementary material hosted on biorxiv doesn not seem to include these scripts/parameters.

I'm also curious: on what basis were were parameter combinations judged to be more optimal than another? I like that you explored how these parameters impacted phylogenetic inference, but it's unclear what the actual optimality criteria were.

Given that you are using transcriptome assemblies and presumaby amino acid sequences since you're using the "protein" method, it is important that some degree of protein redundancy reduction takes place prior to analysis with BUSCO (i.e. using CDHIT or retaining only the longest protein per assembled gene). Use of transdecoders '--single_best_orf' will still retain alternative isoforms which will lead the count of "duplicates" by busco to be inflated.

I think a graphical depiction of this may be useful. As written, it's unclear how you are encoding orthogroups within homogroups.

Based on your description below, it appears to me that the only difference is in the implementation of MCL clustering to sequence similarity, with orthofinder2 used to infer orthogroups (i.e. normalizing similarity scores to account for sequence length), and homomcl for homogroups). Your statement that orthogroups are nested within protein families doesn't appear to be consistent with this? Am I mistaken in my interpretation?

Could this not be said about gene content evolution as well? Patterns and dynamics of gene family gain/loss is certainly expected to be a heterogeneous process - does the reversible binary substitution model allow for site heterogeny?

I am delighted to see such a thorough exploration of the sensitivity of phylogenomic inference to these parameters, as defaults are typically used without consideration but can, as you say, have profound influence on downstream outcomes.

This is the same type of gene-content encoding used by Ryan et al. 2013 as their second dataset, which coincidentally supported the hypothesis of Ctenophores as sister.

I'm curious as to why an alternative recoding of these gene families was not considered - gene copy number. From the analyses conducted here to infer orthogroups using OrthoFinder, this encoding would be fairly straightforward to accomplish, and could perhaps be even more information rich than presence absence. Additionally this would be an entirely distinct way of encoding these data as compared to the approach employed by Ryan et al., (2013).

Additionally, there exist now several methods that are capable of inferring species trees from multi-copy gene family trees, which are by design more robust to, and even improved by the existence of paralogs (e.g. Asteroid -https://doi.org/10.1093/bioinformatics/btac832, SpeciesRax - https://doi.org/10.1093/molbev/msab365, ASTRAL-Pro 2 - https://doi.org/10.1093/bioinformatics/btac620).

Unclear a priori why this might not always be expected to be the case, even when encoding data in other ways. Fundamentally, these are still evolutionary models.

This is an impressive paper containing a significant amount of work, and I love the amount of effort focused on exploring how sensitive the results are to methodological approaches and parameter specifications. In particular, I'm delighted to see the exploration of how the choice of inflation parameter impacts downstream phylogenetic inference - given that this parameter has profound influence on all downstream analyses, I think all studies using these types of approaches should be similarly cautious/considerate.

In general, I feel that the choice to encode 'gene content' as gene presence makes distinguishing these results from past efforts, particularly that of Pett at al., 2019, more challenging. That is, gene content is encoded in the same way across the two studies - and both the present study and that of Pett et al., come to the same conclusion that the Porifera-sister hypothesis is most strongly supported. Together, it makes these new phylogenetic analyses of gene content less compelling as corroborating evidence.

I feel that an alternative encoding of gene content - as gene copy number as an ordered-discrete character e.g. (0,1,2,3, etc) - would have been valuable to explore. This type of encoding would be more information rich, and would be inherently distinct from past efforts. Alternatively, the authors could have used a growing number of methods that are capable of inferring species trees from multi-copy gene families (see later comment). I would hope to see some discussion of these alternatives, though i suspect implementation of either would be non-trivial.

Could you clarify what you mean by this? It's difficult to intuit what this discordance looks like without the supplement being available. Regardless, it's nice to see that the discordance between the two topologies seems modest -

It might be really interesting to see the distribution of D-statistics pooled into both within vs among species and ecomorph comparisons.

I say this because it could be really informative to compare/rank the evidence for hybridization across levels - within/among species, as compared to within/among ecomorphs.

Oh, I understand now! It might be worth adding this explanation in the methods? I would also suggest including a plot of LogLik for each value of K - it's difficult to evaluate how strong of a difference there is in model supports without it.

However, it is somewhat surprising that that these data would suggest there are only 2-3 clusters, when there are 15 described species - this seems like quite the discrepancy.

If i'm understanding the description of the UCE tree correctly, would it be correct to say that there is one additional mitonuclear discordance when considering the UCE tree? (The placement of T. anuenue)

Which nuclear tree? I see below that there is some discordance in topologies inferred using the two nuclear datasets (SNPs with fastme, UCEs with IQ-Tree), so it would be worth specifying here.

Just a suggestion, but it might simplify things for you if you had a section prior to any use of ANGSD describing the commonly used parameter settings across these analyses?

Could you elaborate on why you chose to use this tree versus the UCE likelihood tree? Are results from Dsuite robust to the choice of species tree used?

Was this quantified somehow? I might have missed it, but I don't see this described in the methods.

Could you clarify why 15?

This is a very specific number to use as a filter. Was this determined based on patterns of LD-decay inferred from the steps above? I might provide some additional explanation as to why this was used as the LD-filter.

I have a couple thoughts about this - is there a reason why a combination of fastme and RAXML was used as opposed to full likelihood inference using RAXML or IQ-TREE?

I ask because both RAxML and IQ-tree include methods (that i would suggest using here) for acquisition bias correction when using SNPs that could improve inference of both topology and branch-lengths.

Another entirely different strategy to infer a species tree in a computationally efficient way using SNP data under the multispecies coalescent would be to use SVquartets (Chifman and Kubatko 2014). This potentially well suited for your dataset as it can leverage a large number of SNPs effectively, and even robustly handles linked sites. Caveat is that it cannot infer branch-lengths, but this might not be a concern for your use case.