While this may be true, I'm not sure what you have shown can disambiguate the limits of "traditional" or "standard phylogenetic approaches" from the inherent limits of using RNA-seq as a biological readout and/or the complexities of this type of study design. How confident are you in the fidelity of the transcriptomic profiles or your power to detect joint species effects, given the number of pangenomes and phylogroups? Turning to the protein embedding space is an intriguing idea with some clear utility, but I wonder if a deeper exploration of the relationship between phylogenetic and measurement limits might help build a stronger case for which types of alternative methods are needed.

What are your thoughts on the relationship between µPULLI-I features and µPULLI-DL representations? It seems like a deeper reconciliation may be possible and informative.

Both methods are inferring some latent biological structure that they leverage for prediction. µPULLI-DL is doing this via its pre-trained features; µPULLI-I via pre-defined measurement sets. The ultimate goal, I assume, is to id the minimal dimensions characterizing the physiological state of cells. This need not map directly onto the µPULLI-DL or µPULLI-I features but they are likely useful guides along the path.

Have you considered using tools from the mechanistic interpretability world to probe whether or not the µPULLI-DL features map onto the µPULLI-I ones?

How much variance in the 70-dimensions do the first 2 PCs account for?

The discrete structure of the PC space is interesting. I wonder how the space would change if you were to train the model on non-instantaneous windows. I would hypothesize a smoother space but, maybe paradoxically, the possibility of greater classification accuracy given the autocorrelation of the velocity components.

Also, the 70-dimensional latent space post-optimization is intriguing. Any idea why that size is optimal with respect to the input features?

It seems likely the foundation model embedding geometries could vary across multiple axes: distribution of (non)linearity, local vs. global dimensionality, etc. I wonder if it might be worth characterizing these as a prior for intertepreting the proposed metrics. For example, won't the Euclidean distance of a gene in a locally nonlinear region of an embedding mean something different than a gene in a linear portion?

These results are intriguing but largely observational and/or rely on pre-determined cutoffs. Given your tree and your set of features, shouldn't more quantitative resulsts be available? e.g., are there detectable signals of feature covariation across evolutionary time? hierarchical/directional relationships between directons and genomic features?

I don't think you have yet shown enough to make this claim. The paper would greatly benefit from the use of null models. Given that all the data used in the paper share informational structure (i.e., they're way detectable in some way as "songs" to humans), are the relationships you identify actually surprising? It's very hard to gauge this without a null, either empirical or simulated, to compare to.

Similarly, it's impossible to know whether the comparisons are fair without a deeper examination of how parameter choice (e.g., window sampling size) may differentially affect the CES estimation across song types. Even if you find significant similarity between human/animal song times compared to a null, how will you know if this isn't a product of bias in your CES estimation?

Why give developmental pathway differences such a privileged role here? Developmental expression variation certainly contributes to phenotypic variation, but the same argument could be made for other, potentially more relevant, levels of biological organization. For example, the diversity of vertebrate neural architectures/functions can be found among ray-finned fish. Neural phenotypes also rely on "pathway expression," but, at least in my estimation, are much closer to the explanatory level needed to understand constraints on collective behavior. Fish behavior invokes a number of neural structural and computational innovations that have been cataloged. Transcriptomic + genomic data are easier to generate, so, of course, we have more readymade hypotheses built on them. However, it feels somewhat unnecessary to restrict explanatory hypotheses to these when others are potentially available.

I'm having trouble parsing the distinction here. What would a direct gene–phenotype relationship look like? A phenotype that varies due to a single, highly penetrant gene? Couldn't "distributed integration" also encompass the effects of such a "direct gene–phenotype relationship" in the form of master regulators or molecular information bottlenecks? Presumably, a non-functional master regulator gene would impact distributed integration and represent a direct relationship.

The relationship between what is measured, how it's measured, and evolutionary diversification is always important and interesting. It would be useful to know how the previous aggregation metrics covary with yours. Are they simply more and less nuanced measurements of the same trait? Measuring fundamentally different things? Are they differentially sensitive to certain strains over others? Do the previous metrics show any association with genetic/evolutionary factors? I wonder if this might help refine the true phenotypic locus of selection.

Ancestral sequence reconstruction confidence famously decrease with evolutionary distance. Fitting with this, the posterior distributions in Dishman et al. are less confident for older ancestral proteins.

ASR can also be influenced by the tree and method used. Dishman et al. used PAML on a tree containing 457 chemokine sequences.

Have you looked into the effects of using a different tree and/or ASR method on the reconstructed sequences and their posterior distributions? One worry is that different approaches might product ancestral sequences more/less similar to those AF has seen which, in worse case scenarios, could propagate bias through downstream analyses.

Increasing the ancestral sequence # using different methods might help give intuition around this....

It would be useful to include a bit about the criteria used during expert evaluation

It looks like there might still be some oscillatory patterns as a function of run order in panel 2a. It could be useful to plot a spline for quick visual inspection. Also, I wonder if a window- rather than a point-based correction might better remove this confound?

How worried are you about data leakage influencing classification performance? Does classification accuracy vary with the sequence ID% % used?

How were the clade sizes estimated? I would guess that the accuracy of size estimates will vary quite a bit. For example, there may be accurate estimates for some Eukaryotic lineages but for less-resolved Bacteria we may be way off the mark. We don't know what we don't know in many cases and there could be a good amount of latent variation unaccounted for (especially among Bacteria).

This seems important since the measurement of discrete events (e.g., gains) may be influenced by taxon representation (i.e., larger taxa have more opportunities for detecting gains).

It's hard to parse how much phylogenetic diversity BaseData actually adds (and how it's distributed) with just these plots as evidence. Quantification would be useful.

What proportion of phylogenetic diversity in GTDB/OMG is captured? How much new phylogenetic diversity is contributed? How many novel clades without GTDB/OMG ancestry are added? Crucially, how are these patterns distributed within BaseData? Are phylogenetic diversity gains evenly distributed over proteins/genes?

It's also notable that Phyla just barely outperforms Hamming Distance for the TreeBase dataset but substantially outperforms it with TreeFam. Might this be related to the inherent differences between reconstructing species vs. gene family trees (species trees are often estimated over sets of gene trees)? It could be worthwhile to consider species and gene tree reconstruction as different classes of tasks, one likely much harder than the other.

Since Robinson-Foulds has some well-known limitations, it would be useful to compare it to other distance metrics (e.g. quartet distance).

It might be worth noting that another limitation of tree reconstruction is that it's inferential and probabilistic; actual "ground truth" relationships are never available. Trees are estimates of relationships given a set of parameters and assumptions. Different methods will estimate relationships differently.

What are your thoughts on the possibility of homology-based data leakage (e.g. https://www.biorxiv.org/content/10.1101/2025.01.22.634321v1.full)?

Given the hierarchical relationships generated by evolution, many genes will share some degree of relatedness both among the "phylogenetically diverse species" and with the held-out species. For example, a deeply conserved metazoan gene with little sequence diversity will be a nonindependent data point; it's sequence will be pseudoreplicated (possibly its expression too), allowing information leakage between pretraining and validation sets.

It seems likely that data leakage will happen to varying degrees based on the conservation of each gene. If this is happening to a substantial degree (i.e., the model is learning high copy number/homologous/pseudoreplicated sequences best), then it wouldn't be surprising that performance would scale with evolutionary distance from humans; the amount of shared homology would predict model performance. It also wouldn't be surprising that TF-Metazoa outperforms other models; by having more pseudoreplicated sequences, it provides more opportunities for data leakage and, thus, overfitting.

Is there a convincing way to show that this isn't the case?

How were these species selected? Are there reasons to believe that they're representative of general patterns? Are there any species that don't match known biological constraints? Understanding what the underlying distribution is, of which these species are samples, would be very informative.

In general, it would be helpful to the reader throughout the manuscript if justification were included when presenting results on subsets of species.

Taxonomic biases of the training data likely also play an important role here. The data sources aren't equal in how they've sampled the protein universe. This is especially apparent when comparing the structure and sequence databases. For example, certain taxa (e.g., humans) are overrepresented in the PDB, while others dominate UniRef.

It's not hard to imagine how the distribution differences in the t-SNE might reflect this, especially given the strong overlap of sequence-based methods with the UniRef samples. Do you know if the same is true for the structure-based methods? If you visualized where, say, PDB proteins are, would there be strong overlap?

Any ideas on how to disentangle approach from the taxonomic makeup of training data?

Did you perform any comparisons other than Passeriformes vs. not?

An examination of the 3 parameters (speciation, extinction, diversification) as a function of all taxonomic comparisons seems like a useful and potentially more agnostic analysis. More generally, I wonder about the extent to which taxonomic coarseness might influence sensitivity to detecting "cognitive buffer" over "behavioral drive." Might there be smaller clades of large-brained species that defy the macro-level extinction trends?

ESM2 performance can be sensitive to the makeup of training data used (e.g. https://www.biorxiv.org/content/10.1101/2024.03.07.584001v1.abstract). Specifically, class biases in training data can be recapitulated in generated sequences.

Given that AUC-ROC varies as a function of compartment type (Fig 1D) and the compartments themselves are associated with diverse input sequence numbers (Fig 1B), I wonder if you examined possible biases in ProtGPS's behavior? Does ProtGPS more readily generate sequences that are suited for certain compartments than others? Is this explainable by the statistical distribution of the training data?

There's some evidence of sex-specific cell type heterogeneity in organs (e.g. https://pmc.ncbi.nlm.nih.gov/articles/PMC10210449/; https://pmc.ncbi.nlm.nih.gov/articles/PMC7615307/#S1). It seems possible that consistent sex-specific organ heterogeneity might be another explanation for the patterns you see and, if present, could change interpretations/conclusions. E.g., sex-biased differences could arise from cell number variation rather than intrinsic transcriptional differences. How much of a concern is that here?

Does mapping percentage or gene count vary with phylogeny and/or ecology? Put differently, is there any reason to worry that technical variation here might influence your sensitivity for detecting cell type abundance, especially given the low number of replicates per species?

Can you expand on why you chose blastp? There are a number of other (likely more sensitive) alignment methods. Given that many of the analyses in this manuscript rely on specific assumptions with respect to evolutionary age, it seems that identifying the most accurate approach possible would be useful.

Also, why use this specific e-value threshold for all proteins? Proteins often vary in e-value distributions due to differences in sequence length/composition, evolutionary history, etc. Methods that account for this (e.g. OrthoFinder) might be worth exploring.

Is C. pacifica's capacity to live in higher salinity environments accompanied by variation in their motility patterns with respect to non-extremophiles? Given that the Reynold's number varies with salinity, it might be enlightening to measure C. pacifica's speed distribution at different salinity concentrations. I wonder if these experiments might uncover more interesting axes of diversity within C. pacifica that differentiate them from species like C. reinhardtii.

In the manuscript, "fitness" generally refers to landscapes learned by pLMs. However, at other times, it is used to describe the actual landscapes traversed by evolution (via processes like natural selection). Given the limitations of pLMs - including those you cover in the introduction - it feels dangerous to conflate these two. It is far from established that language models are able to infer the true structure of evolutionary processes, much less model the complex activities of natural selection.

This feels important to note since the discontinuity between fitness and trait distributions has been recognized for a long time (e.g. Fisher 1930). Many factors contribute to this relationship, both at the individual gene/protein level and at the level of genetic interactions. It is likely that variation in relationships between pLM fitness/activity will also be affected by multiple such factors (as evidenced by the differences observed even here across the 5 proteins of focus). It is also likely that these will at least somewhat differ from the factors influencing empirical fitness landscapes. Delineating these differences clearly seems to be useful for future model development/refinement.

I find it hard to interpret the importance of these results without more context.

For example, how surprising is it to see this enrichment given the initial n of natural and generated proteins?

How might decisions with respect to tree construction effect the branch length distribution? It seems possible that you would get different a different outcome if you varied the mmseqs parameters or implemented different criteria for choosing representative proteins.

Furthermore - though novel phylogenetic groups are distributed throughout the tree - it would be interesting to know if the overall distribution across clades is predicted by the abundance of natural proteins across the tree. I.e. do clades with more natural proteins in the training data tend to produce more generated proteins?

What was the rationale for using a tree from timetree.org as opposed inferring one from the gene families?

I imagine that a comparing the effects of using timetree vs. an inferred tree on CSUBST outputs would be enlightening. Such a comparison could be an empirical way to assess the effects of topological error in this data set (and would be a nice complement to some of the analyses in Fukushima and Pollock).

It's great that even with the use of fixed rates you see a substantial increase in fraction of bias explained. Since mutation rates obviously do vary, I wonder just how much better you might do using a model that doesn't explicitly fix them...

What does it look like if you just use the branch lengths from the phylogeny to do this weighting? I would guess you get at least some increase in the Spearman correlations and it's a straightforward approach.

It is hard not to notice that the distributions for certain states (e.g. Meerkat vigilant/resting state) are noisier than others. It would be interesting to see a comparison of the variance of these distributions as a function of species and/or state to see if the claim in this sentence is statistically supported.

How much history/ecological data are there available for these species? It could be interesting to pair the phylogenetic patterns with other trait data to explicitly test different evolutionary hypotheses. e.g. is there a relationship with prey type? substrate? depth? biotic diversity?

Do these behaviors vary at all as a function of what prey are used? I'm guessing you tested squid and crabs with P. carolinus as you did with P. evolans?

Presumably motile (squid/crabs) prey would give off a different set of cues that less/non-motile prey (mussels)? Specifically, I wonder if there is a tradeoff between chemo- and mechanosensation that is dependent on the amount of movement? Examining this relationship could be a potential route into the neural computations underlying digging behavior...

How sparse is this matrix? Given an average of ~1,500 UMIs and ~800 cell-cycle genes, I'm assuming the distribution of expression for the cell-cycle genes is quite distributed/uneven across the cells?

If very uneven, I wonder if some of the cell cycle designations might be driven by sparsity as opposed to canonical expression signatures associated with each stage? One way to parse this out might be to look at the PC loadings using as input to clustering/UMAP/etc. Do any show signatures of extreme sparsity (e.g. binary expression only one or several genes)?

More broadly, it might be helpful to report the average # of cell-cycle genes detected in each cell.

I wonder if calculating the autocorrelation of coda durations might be a nice complementary measure here. Autocorrelation could give you a sense of the time scale over which the rubatos decay and, seemingly, might also provide a sense for the timescale of longer-term trends.

Similarly, I wonder if cross-correlation might be useful for comparing the information quantity shared with interlocutors? The correlation value would be interesting, in addition to any patterns of temporal lag between codas. It might be a comprehensive metric for comparing the similarities of codas over time (as opposed to just looking at overlapping codas).

It is hard to tell from the text if HDBSCAN is run on the behavioral parameters or on the UMAP output. If the latter, then I would take extreme caution in thinking about the generalizability of the method given the numerous issues with clustering on nonlinear manifolds. Either way, it would also be helpful to report more information on what the specific morphological/motility parameters are and any normalizations/manipulations that were done on them prior to UMAP and clustering.

Also, any justification for choosing of UMAP and HDBSCAN would be useful.

This might be a slightly confusing way to refer to UMAP dimensions (is it accepted that a UMAP dimension = a single 'UMAP'?)

Does donor identify have any effect here? Do the donors differ at all in their speed/area distributions and effect sizes? This would be useful to know here and for many other analyses presented in the manuscript. More broadly, it is a little hard to assess the generalizability of the behavioral results presented here (including the cellPLATO analyses) without knowing more about the influence of experimental variables like this.

These species are ecologically diverse (e.g. benthic vs. riverine) and likely possess corresponding sensory differences. Given this, it seems possible that their prey capture behaviors may vary as a function of sensory environment. For example, benthic species may display different repertoires in dark conditions.

Have you tested the effect of varying the sensory environment on prey capture behaviors? Is there intra-specific variation? Are species-specific behaviors invariant? Whatever the outcome, these experiments would help refine the picture of how these behaviors evolved and could lead to more specific sensorineural hypotheses.

I wonder if it might be worth including a brief comment on the identity of these genes and/or their potential relationships with echolocation? Do they seem to be sensible functional hits? Seeing as the echolocation score appears to work quite well it would be interesting to known a bit more about any molecular context for these predictive loci.

phenotypes

I'm wondering about how varying these criteria would effect the number of/which genes were detected.

For criteria 1, what was the rationale for choosing agreement between >2 algorithms? From Fig 1C, it's hard to tell if there is a relationship between %homology and #of agreeing algorithms. What benefit do you get from using this cutoff? What are the tradeoffs? It might be helpful to include a figure similar to 1C, but including the full set of genes before filtering and to walkthrough the outcomes of different cutoffs.

Similarly, for criteria 2, what type of/how many hits do you get if you don't select for 'neuro' or 'musc'? Is there any chance that, though you are using a behavioral readout, genes not annotated 'neuro'/'musc' might still contribute to a behavioral phenotype (e.g. via pleiotropy/epistasis)? Would be useful to include a statement of your thinking on this!

Seems like it would be worth including a direct comparison of Brownian motion to other evolutionary models. The computational overhead shouldn't be very high and, if the comparison supports the use of Brownian motion, it could be a more compelling argument than this.

I wonder about the effect of scRNA-seq methodology on downstream results here. How do droplet-based approaches (like that used for van Zyl et al.) compare to others (e.g. Smart-seq2) when generating cell type trees? There can substantial differences in the # of genes detected by these methods, with droplet-based approaches often generating datasets with less genes. Does this affect the estimation of rank and/or the outputs of the PCA you use for evolutionary modeling? It seems like this would be an important issue to solve since droplet-based methods are essentially downsampling informative data in a nonrandom way that may bias evolutionary inference.

TLDR: are cell tree topologies consistent independent of sequencing methodologies?

Is there a reason for the coarse sampling at 0.001 Hz? Mechanical constraints of the XY plotting robot? Data size constraints? Obviously faster sampling would open up locomotion/behavior as a read out of other possibly interesting, orthogonal phenotypes (with their own developmental modes). Given the video data you are already collecting, it seems like if faster sampling is possible this would be a relatively straightforward - and informative - set of phenotypes to add in?

Is there a reason for the coarse sampling at 0.001 Hz? Mechanical constraints of the XY plotting robot? Data size constraints? Obviously faster sampling would open up locomotion/behavior as a read out of other possibly interesting, orthogonal phenotypes (with their own developmental modes). Given the video data you are already collecting, it seems like if faster sampling is possible this would be a relatively straightforward - and informative - set of phenotypes to add in?

What about intra-annotator variability? Seemingly this could also be an important contributor to inter-annotator variation. Might it make sense to compare multiple annotations from a single annotator and use the average as the basis for the ethograph generation?

How do these comparisons look for other (potentialyl more 'complex') behaviors? Presumably, turning should be among the more straightforward behaviors for a human to recognize. Do the patterns of inter-annotator agreement change with other behaviors (e.g. grooming) and, if so, would accounting for this increase/descrease performance of the neural network? This is a general risk when using human-based annotations for behavioral classification and it seems to me not easily solved by focussing on a single behavior.

I wonder if it might be useful to do some analyses of the specific spatial orientation of nuclei across the centrifugation experiments. While it is sensible that density may be the primary signal driving cellularization, it is also interesting to consider that there may be higher order relationships between cell distribution or spatial organization that are predictive of the different outcomes (i.e Flip, lysis, irregular invaginations) since centrifugation is a relatively forceful and disruptive approach. Spatial relationships of nuclei could theoretically be extracted by segmentation/registration and performing some basic statistical comparisons to uncover the relationship between the images (e.g. PCA).