Thank you for going through the effort of putting together this properly benchmarked analysis, the lack of a purely additive model in the original work is a significant omission. To give credit to the original authors, it seems the biggest success here is their ensemble PLM model, I suppose some claim of 'synergistic epistasis' could be made here, as the mutants the ensemble proposed do seem to be genuinely beneficial. However, it is obvious that this does not extend to the MLP trained by the authors where the claims of epistasis being captured are made.

One thing I personally have always wanted to see in these sort of data/experiments is some flavour of variance partitioning as in quantitative genetics. I.e. How much variation does 2nd order epistasis explain in these systems? Consequently what even is the expected improvement we can hope if we train a model that can successfully capture say 2nd order epistasis? It seems from these results it would probably be quite modest (at least in the context of these data). I don't think a proper analysis can be done here due to the experimental design; each mutant combination is measured only once, so epistasis cannot be distinguished from experimental noise. This paper from the Thornton lab has done an analysis along these lines, and the results seem to suggest additivity can get you a lot of the way there...

This is super cool work leveraging conservation across deep evolutionary time scales. The core thesis is that the PPI that are old (and detectable) enough across these time scales should be fruitful for genotype-phenotype mapping in contexts like human disease, makes overall sense. The fact that you can point to specific examples where this works is also quite remarkable. However, I do wonder if there's potential interesting followup to more directly test the hypothesis of ancient homology -> relevance to disease. The approach outlined here may be just as useful on shorter (but still deep timescales) e.g. across vertebrates, where homology/conservation might still be informative, but there might be more ability to pick up signal on these sorts of relationships. It would be extremely interesting to see does e.g. Fig 6A changes as you alter the timescale of divergence in the selected species. Going shallower will likely come with the potential drawback of adding more noise in the analyses, but this is a relationship that is worth explicitly interrogating.

It would be good to have more detail on how this shared latent representation is achieved. Yang et al describe using a discriminative training objective to achieve this which is not mentioned here. Another even simpler approach is to simply train a linear layer between a frozen input encoder and output decoder.

One note on GWAS modelling which the omnigenic model relates to, is that such analyses require care in dealing with covariance amongst inputs (i.e.in the case of GWAS linkage disequilibrium and population structure). Otherwise polygenicity can be conflated with allelic covariance across many sites. Similarly, it would be useful to get a sense of the covariance structure for both the input and output data in this study. The fact that a linear autoencoder with a bottlneeck of 150 latents has good reconstruction error suggests that collinearity must be fairly extensive. So for example seeing how 'polygenic' your predictive signals are in expression PCA space is might be quite informative. If prediction comes from a few PC's vs many tells you something about how correlation in expression state is distributing learned signal.

I really enjoyed reading this preprint, I think similar lessons on what exactly CV evaluates depending on your strategy have been (re)learned several times across deep learning biology.

typo

typo

There's also other secondary costs to aggressive recurrent GS such as drift in traits that aren't amenable/feasible to build GP models for.

I understand the rational for why the current study is setup to ignore CV0/CV2 scenarios, however I would argue that CV2 is likely one of the more promising study designs for genomic prediction at least for commercial breeding programs. By leveraging sparse testing + GP, one can test material across a broader environmental footprint for a fixed budget. The GP based breeding values from such a study design as particularly valuable as they can allow more relevant selection decisions to be made, as commercial breeding programs are often aiming for broad acreage products as the end goal.

It might be valuable to check which conformation (if not both) were included in the original model training datasets.

It's interesting to see how the key residues in these attention maps interact globally with the total sequence. This feels somewhat distinct from the results of Zhang et al. on the categorical Jacobian which picks up strong pairwise patterns between amino acids (predicting the contact map of a folded sequence). I wonder if this pattern is a unique feature of these fold-switching proteins or a general phenomenon in Alphafold.

Very cool approach!

Might be another typo in the figure here "E2P"

This is probably the second most informative model comparisons in the manuscript in terms of validating the BINN framework after the permuted gene/eQTL BINN models. It would be helpful to see a statistical test of predictive power between the random vs bio-informed ridge regression models.

I think it might be useful to add an asterisk to this statement. In principle by tracking the cascade of G->B->P one might construct a more powerful model through biologically informed sparsity. However, it is crucial to keep in mind that while expression variation is partly driven by genotype, in a field setting realistically a significant fraction of expression might be attributable to environmental effects. (This tracks with the performance of expression only models in Sup fig 3/4). E.g. a stressed plant likely has a unique transcriptomic profile, and an easier to predict silking phenotype. However, these relationships won't translate to the G->B->P model since transcript abundance itself is not used during model inference. This effect is what likely explain the difference in performance between the expression RR models in sup fig3/4 and the G2B2P model.

The results from Fig 2 in the main text do not seem to suggest that the B2P model is substantially more powerful than the G2P model which might suggest that this could be an issue in the maize dataset. On the other hand the Fig 3/4 in the supplement do show (the expected) pattern that expression is the most powerful predictor of phenotype. I have struggled to reconcile the results from these two figures while reading the manuscript.

This doesn't match the "Anthesis MI" label in the figure, might be a typo.

One more dimension of robustness that I think could be useful to explore is sources of error in data (e.g. genotype/phase etc.). Actually adding some small amount of noise to genotypes during training could make cxt quite robust to real dataset errors giving it a further performance edge over alternatives like Singer.

In principle do you reckon it would be possible to use a random masking approach (like in models like ESM2) for this problem? Currently one (entire) side of a region informs the models prediction on the focal window, while in principle the most informative regions are both to the immediate left/right of the focal window. Random masking as a strategy could allow the model to leverage this information bidirectionally, but technically could be more challenging.

We greatly enjoyed reading this preprint for our internal journal club. This seems like a very principled and useful application of the transformer architecture in biology.

It would be very interesting how this effect scales with model size. I'd be curious to see if particularly the larger ESM models cause more general latents to break down into more family specific groupings. This would mesh with some of the evidence of family specific overfitting that keeps popping up in the larger sized models.

This way of backing out recombination rate from Ner makes some strict assumptions, particularly that pi reflect mutation-drift equilibrium, and that mu is invariant across the genome. This is unlikely to be true in most datasets. It would be very helpful if more region specific comparisons could be made to pedigree based estimates to evaluate how well this approximation works in this dataset.

Again I think the influence of phylogenetic signal is important to consider here. A phylogenetically corrected regression (e.g. phylogenetic generalized least squares) would take into account the non-independence of species observations and provide a more accurate view of the data.

It would be nice to get an idea of how stable this porting over of the estimated recombination map is across phylogenetic distance. There is considerable variation in the correlation among species pairs which would be good to understand. Also would be valuable to show the reader how the phylogenetic distance of these comparisons relates to the distances used for comparisons of the main analysis.

The message of this figure get's a little muddled. One would expect that as as Ne increases, that the proportion of mutations experiencing a larger (more negative) Nes would grow (as in the case of increasing r) but the opposite is shown. Unless the figure implies that fewer segregating mutations will have large Nes? But then this conflicts with the sequence for recombination.

I understand that for computational efficiency/interpretability it is probably best to restrict the Borzoi comparisons to a difference in one focal variant. But given the size of these windows there will almost certainly be many variants that differ across sequences in the study population. I'm curious if you have experimented with using actual UKBB haplotypes as input for a focal position and tested if this introduces meaningful variance in the predictions for said focal variant. Could be valuable to assess how much artificially restricting sequence space differences to one position affects model predictions.

The figure reports an average but my gut feeling is that this threshold should maybe be protein/protein family specific. I imagine the overall shape of the phylogeny/distribution of branch lengths around a focal protein will influence how much predictive gain LFB provides. For example, it might make sense to set this threshold higher for a protein with lots of intermediate divergence homologs, vs one that has few. An explicit analysis of what features of a protein family's phylogeny favour differing thresholds might in itself be a very useful analysis for guiding the application of LFB/LFB-like methods to PLM improvement.

One alternative/complementary 'modify the model' strategy that might be useful to compare to this method is protein family specific PLM fine-tuning. One could test how fine-tuning on a much narrower region of protein space affects a PLM's ability to soak up phylogenetic signal by testing if a fine-tuned model is similarly improved with LFB in zero-shot fitness prediction tasks.

With an emphasis on model generalizability, the most interesting signal one could observe is that embeddings cluster not by species (which is driven by strong phylogenetic signal) but rather by other conserved biologically meaningful differences like cell type. An explicit quantification of how much clustering in UMAP (or some other dimensionality reduction method) is explained by species identity vs cell type would be convincing of model generalizability (and the benefit of having multiple species in the training data). At a glance the plots right now suggest the dominant driver in clustering does seem to be species identity but it is hard to tell.

Is there some quantification of this? In the challenging cell types for example the performance seems roughly equivalent between the Metazoa and Exemplar versions of the model. Overall it is hard to see evidence of a benefit of adding diverged species to the training data.

It would have been interesting to see how a simpler say vanilla MLP based approach would stack here to really sell the advantage of attention over other deep learning approaches.

I greatly enjoyed reading this paper. The rigorous and rational approach to testing model performance on simulated data, reasonable model architecture, and smart dataset choice are a much needed advance beyond haphazardly applying deep learning networks to G-P datasets with minimal performance gain.

Did you experiment with how LD based pruning affects model performance? For linear genomic prediction models the relationship between marker number and predictive performance is well characterized, as long as you capture LD structure well, marker number is not very critical. However, this is not characterized well for deep learning models in this context. Epistatic interactions in particular will depend on products of LD between marker/causal QTL's which could cause performance degradation if causal QTL's are not very well tagged.

I was curious why you chose to simulate fourth order epistatic interactions. Statistically one expects higher order epistatic interactions to contribute progressively less to genetic variance, so most studies tend to focus on pairwise epistasis.

It doesn't seem that Fig S5 shows how layer 26 was selected. It would be interesting to at least get a short description in the methods of how this layer was chosen. Other work on mechanistic interpretability in protein language models has shows that different types of features can be learned in different layers of the model.

As an evolutionary geneticist to me the most interesting benchmark here are the PhyloP scores. When I see models like EVO2 my concern is always that they are able to effectively memorise phylogenetic conservation. This is totally valid from a biological standpoint however, this can be done with a far simpler phylogenetically explicit method like PhyloP, GERP etc. What is far more exciting is the possibility that a flexible, large model like EVO2 could pick up on non-linear (e.g epistatic) patterns which is something PhyloP type methods are blind to. That PhyloP is very competitive in all these tasks I think is quite telling that for the most part the power of all these models comes from identifying conservation rather than more general 'biological rues'. However that in some instances PhyloP can be improved upon is very exciting nonetheless, in my opinion this is the golden benchmark to be trying to beat.

Aside from NT, these alternative metrics of lncRNA essentiality seem over simplistic compared to a model as complex as EVO2. Are there no other alternative models for lncRNA essentiality? Maybe a tweak of sequence conservation methods could work here too.

The results certainly seem to suggest these native populations might be bottlenecked too. Is there any indication on how central these sampling locations are to the species native range? Is it possible that the range edge was sampled?

Site-frequency spectrum plots per population+mutation would quantitatively demonstrate these patterns without the need to arbitrarily bin allele frequency.

This correlation is based on two autocorrelated measures (as theta pi synonymous is measured in both the X and Y axis), so it should be interpreted with caution.

I wonder how much the fact that DVV is a crop pest might influence the results for this species. It would be easy for me to imagine that most DVV populations (native and invasive) have experienced agriculture related bottlenecks and/or population expansions. Pests like corn rootworm have repeatedly adapted to the use of pesticides/GM crops a process which often involved a bottleneck (followed by expansion) and may cause similar effects on the evolution of load in native/invasive populations. Data on the population ecology or local agricultural practices (and history of pest load) may be helpful in figuring if the selective landscape of these populations could have such effects

Very cool result!

It would be helpful to have colour labelled legends in the figures.

As you point out, it seems some of the strongest evidence for ECG is that polymorphism is not just shared, but shared polymorphisms are linked. One way you could statistically quantify this is by running a tool scanning for evidence of identity by descent (IBD), or use a tree sequence approach, treating each gene from each accession as an individual genome (like in the multiple sequence alignment you construct). This isn't strictly what IBD tools are for, but it should provide a good proxy given that A. thaliana has low polymorphism and high linkage disequilibrium. Relatedly it would help if intervals for putative EGC could be filled in, not just the limits marked as in Fig S6. This would make it easier to see what the length of EGC tracts is.

This seems to me as a fairly strong statement that doesn't quite line up with the goals/results of this study. Technically what this study shows is how ECG can contribute to standing genetic variation in a population. The specific paradox of standing genetic variation in phenotypes however, looks to reconcile why there is more variation than expected under mutation selection balance (MSB). MSB in practice is agnostic to the type of mutations, when/where they arose, simply their fitness effects. As detailed in the first reference, it is clear that SNPs alone are inconsistent with MSB which is not surprising since they are only a fraction of genetic variants found in populations. However again as detailed in the first reference, quantitative genetics approaches that use mutation accumulation experiments are agnostic to mutation type, and provide a framework for testing if MSB is sufficient to explain standing levels of genetic variation. Paradoxically they often find that MSB is not enough to explain the high variation in phenotypes (due to genetic effects) we observe (see https://doi.org/10.1098/rspb.2018.1864 for an example), implying that forces such as balancing selection must also be working to maintain excess genetic variation. What this study demonstrates is that non-SNPs can contribute to variation (as other studies have demonstrated for transposons, inversions, indels etc.). This alone does not demonstrate the sufficiency of MSB to explain observed genetic variation in phenotypes though.

It might be interesting to see what the scale of hamming distance distribution is in the underlying MSA's for the focal protein families is, vs. at what scales of hamming distance such effects are observed in the simulations. One potential concern could be that couplings/epistasis are estimated from the MSA on one scale of sequence divergence, but the simulations are pushed to much larger scales, in which cases the epistatic interactions inferred from the MSA might no longer be accurate.

I enjoyed reading your paper! This canid dataset offers such a great opportunity for exploring genotype-phenotype mappings, it's great to see how such associations can be teased apart in studies such as this.

Relatedly, it might be interesting to compare effect size estimates across these different data subsets. A large swing in additive effect sizes for markers across populations has implications regarding epistatic interactions the focal locus may be part of, see: https://www.nature.com/articles/nrg3627

It would be interesting to see a simple summary in the text of how much sharing there was in significant markers between the breed average/small/large subsets.

It's hard to tell from the description in this manuscript how non-additivity is captured by this analysis. A quick one-liner on this might be helpful for readers.

Do you have a sense for how much of a difference there is between using an ROH sharing vs a general whole-genome SNP/marker based kinship matrix? My initial reaction was that a more independent whole genome derived matrix that captures structure across the whole genome might be more desirable to capture fine scale population structure/ancestry differences etc. But perhaps using ROH runs themselves since they are the target of the analyses is better.

A small content suggestion. The introduction of this paper spends a considerable amount of time discussing the potential history of dog domestication. However, this background seems only tangentially related to the content of this paper which rather aims to take advantage of population structure in dogs to explore associations between genotypes and phenotypes in the context of domestication. More background on the genomics of domestication might be more relevant.

I'm a little confused by this statement. This would make sense in certain cases (e.g., highly recessive mutations). But SNPeff and GERP scores give no information on recessiveness. If these are just any deleterious variants, then the expectation should be the opposite, that high Ne pops should purge load easier.

Does this seem plausible in the case of kelp? If bottlenecking has been recent, the effect on Ne will be instant, but it will take for the signal of differences in purging to build up. Have strong kelp declines been more recent (>50 years)?

Seems in bull kelp lots of gene flow occurs across the southern tip of Vancouver Island. On the other hand, the northern tip seems to represent a barrier (judging by the clusters in Fig1). Are there any hydrological/oceanographic reasons to expect this maybe?

I really liked reading this paper. It's great to see such detailed sampling and interrogation of the pop-gen of these keystone species.

Does this mean you report selfing Ne (rather than typical coalescent Ne) in your analyses? This would be useful to highlight in the main text.

Is there any correlation in Ne/inbreeding/low diversity between the two species of help where they co-occur? This could be a useful indicator for conservation efforts.

While there are no differences between populations, the regression lines for the GERP/SNPeff analyses clearly show that less constrained sites harbor more diversity than more constrained sites, implying that purifying selection is acting on purging deleterious variants in both species. Seems like purging is present in this dataset it is more of a time-scale issue when it comes to detecting it. This makes sense particularly for recessive variants since they will be hidden from selection for a lot of time.

I greatly enjoyed reading this preprint. The dissection of the origin of the de novo contingency locus was very cool.

Do you have any thoughts on why global mutator alleles underpinned evolvability in two populations, and a local mechanism in the other? Seems that increased mutation rates are a common by-product of experimental evolution (e.g. instances in the LTEE). There is a nice paper in yeast that has demonstrated that mutator alleles tend to be favoured in cases where local population size is high (which allows selection to more efficiently act on the beneficial variants they produce). Might be relevant here: Sign of selection on mutation rate modifiers depends on population size, Raynes et al.,2018

What a cool result. This reminds me of the observation in stickleback that the independent evolution of loss of pelvic hindfins tends to target the same locus because of the specific molecular characteristics of that stretch of sequence. This may be of relevance to this study: DNA fragility in the parallel evolution of pelvic reduction in stickleback fish, Xie et al. 2019.

It seem to me that the key component of the experimental selection regime is that individual level and lineage level selection were allowed to act in separate consecutive timesteps, shielding lineage level selection from being swamped out by 'short-sighted' individual level selection. I wonder how likely this scenario is to play out in circumstances such as pathogen evolution where I would imagine that both levels of selection should still be acting concurrently. I agree however that the presence of contingency loci implies that some form of ecological conditions likely exist that allows for this to happen.

Could you constrain this analysis to mutations that are in LE with other de-novo mutations to test this hypothesis?

I greatly enjoyed reading this paper. True experimental estimates of the DFE in MA studies are super valuable and provide a very interesting comparison for pop-gen based DFE methods as pointed out by the authors.

Is there a chance that false negative mutations (i.e. incorrectly unobserved events in the MA lines) could contribute to this result?

Two quick thoughts for further sanity checks. 1) Does this regression look any different for SNPs vs indels? 2) Do the individual mutation specific effects conform to expectations one might have based on the functional annotations available for these mutational events?

The DE approach seems quite powerful especially since it adds a 'benign E' reference line for fitness comparisons. I would love to see how the prediction from this model lines up with true fitness in figure 2 for all lines tested.

It seems from figure 2 and 3, the dominant pattern in the dataset is that of antagonistic interactions (at least in respect to the additive model). This made me wonder two things: 1) Are there are general biological explanations for such a pattern or considerations for why this might be expected? I'm thinking of the GxG equivalent where we know for example that diminishing returns epistasis is a common feature of adaptive populations, and this can be linked to theoretical models of fitness landscapes in the context of Fisher's geometric model etc. 2) Is this the correct biological null model to use? Certainly in the quant-gen world the additive approach would be the go-to starting point, but is this relevant for the context of these fitness estimates? My first gut feeling was that the average null model should be more relevant. Not sure if a pop-gen multiplicative approach is another potential null.

I enjoyed reading this paper and the novel ExExG framing of the study! This is a great dataset, I hope more genomic data can be attached to it in the future enabling even more mutation specific questions to be asked.

It would be useful to get a short description of these models here (aside from the methods) for clarity.