MainA fundamental challenge in biology is understanding the functional significance of genetic variants. Despite tremendous advances in experimental technologies in genomics, determining which variants influence phenotypes or contribute to diseases remains difficult. Evolution provides valuable insights: deleterious mutations tend to be purged by natural selection, although tolerated or advantageous changes may accumulate. Millions of years of evolution therefore provide a genome-wide record of functional constraint for variant interpretation. This idea has deep roots in molecular evolution and comparative sequence analysis6. More recently, large comparative genomics studies have quantified evolutionary constraints across the genomes of hundreds of species7,8.Genomic language models (gLMs) have emerged as a promising approach for extracting evolutionary information directly from raw DNA sequences through self-supervised learning (ref. 1 and references therein). By predicting masked nucleotides from sequence context, gLMs estimate per-site likelihoods that reflect evolutionary constraint without requiring labelled data. These likelihoods have been shown to be effective predictors of genome-wide variant effects1,9. However, gLMs based on standard language modelling frameworks still underperform compared with much simpler classical phylogenetic models on certain variant interpretation tasks—particularly in complex eukaryotic genomes such as in humans and in distal regulatory elements such as enhancers2—even with massive model sizes3,4.A long-standing approach to modelling evolutionary data is to construct multiple sequence alignments (MSAs). By algorithmically aligning homologous loci across biological sequences, MSAs show site-specific patterns of conservation and facilitate the inference of evolutionary preferences for variants at specific positions. Protein MSAs have enabled highly successful models, including AlphaFold10, MSA Transformer11 and EVE12, and have recently experienced renewed interest as scaling single-sequence protein language models has shown diminishing returns13,14,15. Extending this concept from proteins to genomes, whole-genome alignments (WGAs) of complete assemblies from tens to hundreds of species enable genome-wide studies of sequence evolution16,17. Classical parametric phylogenetic models fitted on WGAs—such as PhastCons18 and PhyloP19—have long been staples in variant interpretation for humans and other species. Recent large-scale consortia have further expanded high-quality WGA resources20,21.Therefore, combining the gLM framework with WGA data presents a promising direction, and we recently demonstrated this potential with GPN-MSA22 trained on a WGA of 90 vertebrate species. Here we present GPN-Star (genomic pretrained network with species tree and alignment representations), a general gLM framework that more effectively exploits multispecies WGAs through a new model architecture with phylogeny-aware designs; see Supplementary Information for key differences between GPN-MSA and GPN-Star. This flexible framework enabled us to apply GPN-Star to three WGAs relative to the human genome, spanning vertebrate, mammal and primate evolutionary timescales. GPN-Star achieves state-of-the-art performance across coding and non-coding variant interpretation tasks. The vertebrate model consistently outperforms GPN-MSA trained on the same alignment, whereas the mammal and primate models excel particularly at predicting non-coding variant effects. We further evaluate its use in human genetics by analysing the heritability of complex traits. Importantly, GPN-Star trained on the primate WGA prioritizes variants with an unprecedented amount of heritability enrichment across more than a hundred complex traits23. Also, we find a striking connection between the most informative evolutionary timescale for constraint prediction and the effective polygenicity of traits.To demonstrate the generality of our approach, we apply GPN-Star to five model organisms with minimal tuning and show its effectiveness in assessing variant effects in these species.Alignment- and phylogeny-informed gLMsGPN-Star learns functional constraints by leveraging evolutionary signals embedded in multispecies phylogeny and WGAs. Inspired by classical evolutionary models, it characterizes how genomic positions evolve across species. To achieve this, we introduce a specialized transformer architecture that integrates sequence context and evolutionary information (Fig. 1a).Fig. 1: Overview of GPN-Star.a, A diagram of the GPN-Star model architecture. The input to the model is an arbitrary multispecies WGA window. The target sequences and source sequences are constructed from the alignment window. The source sequences are compressed into clade-level embeddings using attention pooling following the species tree. The target sequences are encoded through a stack of GPN-Star encoder blocks, where the phylogeny-informed cross-attention module integrates information from the source sequences guided by evolutionary distances between species based on the species tree. Finally, a classification layer transforms the encoded embeddings to nucleotide probabilities at each locus in the target sequences, which are then used to compute the training loss or to make predictions on variant impact. H, hidden dimension; K, number of encoder blocks; L, length of the alignment block; NC, number of clades; NS, number of species; NT, number of target species; V, size of the nucleotide vocabulary; w, loss weight; ℙ, nucleotide probability. Q, K and V (and Q', K' and V') in the encoder blocks refer to the components of an attention module (query, key and value). A full description of the model is provided in Methods. b, Application of GPN-Star to the human genome. Three models were trained on vertebrate, mammal and primate alignments, respectively, learning functional constraints at different evolutionary timescales (Methods). c, A summary of performance across the downstream tasks of the vertebrate, mammal and primate GPN-Star models. For ClinVar, COSMIC, OMIM, HGMD and GWAS fine-mapped datasets, the performance metric was AUPRC. For ProteinGym, the metric was the mean Spearman’s ρ across assays. For S-LDSC, the metric was heritability enrichment. The performance metrics were scaled linearly to the range of 0 to 1; in each task 0 was defined as the lowest metric among the three GPN-Star models and PhyloP and PhastCons fitted at the same three evolutionary timescales, and 1 was defined as the highest metric. cLLR, calibrated log-likelihood ratio; Ma, million years ago. Animal icons in a,b created by Berkahicon are from Flaticon (www.flaticon.com).GPN-Star uses an encoder-only architecture trained with a masked language modelling (MLM) objective. Each input consists of a WGA window and its species tree. The alignment is partitioned into target and source sequences and the model predicts masked nucleotides in the targets conditioned on both their sequence context and the evolutionary information in the sources. This is accomplished through a stack of encoder blocks, each comprising a sequence-wise self-attention module to encode intrasequence context, a phylogeny-informed cross-attention module to encode evolutionary context from the source sequences and a feed-forward network to integrate the information. To efficiently represent phylogenetic structure, we group closely related source species according to their evolutionary distances and pool their sequences into clade-level representations. The cross-attention module adaptively weights source sequence contributions using evolutionary distances derived from the species tree using an attention mechanism. The resulting embeddings for the target sequences are then passed to an MLM head for nucleotide probability prediction. Full architectural details are provided in Methods.Another methodological innovation of this work is that it accounts for mutation rate variation in gLMs. Because GPN-Star is trained on genomic sequences observed in nature, its predictions reflect both mutation and selection. To disentangle these signals and isolate selective constraints, we normalized the prediction by the model of each variant against its background mutation effect estimated from high-confidence neutral sites in the genome (Methods). Unless otherwise specified, all analyses in this study use these calibrated scores to capture functional constraints with minimal mutation rate biases.The GPN-Star framework is general and flexible, designed to work with any alignment data from any species, requiring minimal hyperparameter tuning to achieve robust performance (Methods). We first applied it to the human genome, training three separate GPN-Star models using the largest available vertebrate, mammalian and primate WGAs (Fig. 1b; Methods). Throughout this study, we denote these models as GPN-Star (V), GPN-Star (M) and GPN-Star (P), respectively. We have experimented with different model sizes, but we focus here on the largest models (200 million parameters), trained for several days on 8 NVIDIA A100 GPUs (Supplementary Table 1). Existing single-sequence gLMs such as Nucleotide Transformer (trained on 128 A100 GPUs for a month)3 and Evo 2 (trained on over 2,000 H100 GPUs for several months)4 require substantial computational resources to implicitly learn evolutionary constraints, in exchange for flexibility and broader applications. We show that GPN-Star achieves superior constraint prediction with a much smaller resource footprint by leveraging the explicit evolutionary context provided by alignment data.Unlike previous gLMs trained across extremely broad evolutionary distances (for example, spanning prokaryotes to humans), GPN-Star focuses on narrower, more recent phylogenetic distances closely related to humans (Fig. 1b). As we demonstrate below, modelling longer evolutionary histories is not always helpful. Instead, capturing recent evolutionary constraints proves especially advantageous for interpreting certain classes of genetic variants (Fig. 1c).Pathogenic coding variantsA critical challenge in human genetics is understanding the impact of genetic variants on disease susceptibility—an essential step towards improving the diagnosis of genetic disorders, identifying drug targets and realizing the promise of precision medicine. Many computational models have been developed to predict variant pathogenicity in the human genome. Here we systematically evaluate the performance of GPN-Star across a comprehensive set of benchmarks.Our evaluation of coding variants focused on missense variants, the most prevalent class. We first considered classifying pathogenic versus benign variants in ClinVar, a widely used clinical variant database containing expert-curated pathogenicity labels24. We compared our models with major genome-wide variant effect predictors, including classical evolutionary methods (PhyloP and PhastCons) that also use WGA, the ensemble model CADD25, recent single-sequence gLMs (Nucleotide Transformer 2.5B multispecies and Evo 2 40B) and GPN-MSA. GPN-Star (V) achieved the highest area under the precision-recall curve (AUPRC), matching the performance of the protein language model ESM-1b26 while outperforming ESM-2 (ref. 27) and ESM-3-open28 (Fig. 2a). For ease of visualization, only the best-performing GPN-Star model is shown in Fig. 2; results for all three versions are available in Extended Data Fig. 1.Fig. 2: Performance of GPN-Star on human genome-wide variant effect prediction.Classification performance measured by AUPRC is shown in a–c and e–h; bar heights and dot positions indicate the AUPRC computed on the complete benchmark set and the error bars are 95% confidence intervals from 1,000 class-stratified bootstrap resamples. Sample sizes (positive versus negative class) are shown in each panel heading. a, ClinVar pathogenic versus benign missense variants. b, COSMIC high-frequency (>0.1%) versus gnomAD v.3 common missense variants. c, GWAS fine-mapped putatively causal (PIP > 0.9) versus non-causal (PIP < 0.01) missense variants across 65 UK Biobank traits. PIP, posterior inclusion probability. d, Spearman correlations with deep mutational scanning fitness scores for 31 human–protein assays in ProteinGym; each point is one assay. In each box, the centre line is the median, the diamond is the mean, the bounds are the 25th and 75th percentiles and the whiskers are the most extreme values within 1.5 times the interquartile range. e, OMIM pathogenic versus gnomAD v.3 common non-coding variants. f, HGMD pathogenic versus gnomAD v.3 common non-coding variants. g, GWAS fine-mapped putatively causal versus non-causal non-coding variants across 83 UK Biobank traits. h, OMIM pathogenic versus gnomAD v.3 common (y axis) and GWAS fine-mapped causal versus non-causal (x axis) promoter variants. i, Radar plots comparing GPN-Star with PhyloP and PhastCons at the vertebrate, mammal and primate timescales across the seven pathogenicity benchmarks; metrics were scaled linearly to 0–1, defined by the lowest and highest values across the nine models. j, RVAT on 34 UK Biobank quantitative traits (WES; 161,822 unrelated individuals of European ancestry), comparing DeepRVAT with and without the three GPN-Star predictions (V, M and P): number of significant genes (x axis) versus the number replicating previous larger studies (y axis). Each set of three same-coloured points denotes three random initializations.We next considered somatic cancer variants from the COSMIC database29, constructing a benchmark to distinguish missense variants frequently observed in tumours from common missense variants in the general population (gnomAD30). As shown in Fig. 2b, GPN-Star (V) substantially outperformed all competing models, demonstrating strong predictive capability for pathogenicity beyond germline variants.Beyond clinical variants, we also assessed performance on deep mutational scanning (DMS) data, functional assays that measure the fitness effects of all possible missense mutations in a given protein. Across 31 human DMS datasets from ProteinGym31, GPN-Star (V) outperformed all genome-wide models, although it slightly lagged behind the protein-specific models (Fig. 2d).We also compared GPN-Star with two recent missense variant effect predictors, AlphaMissense32 and PrimateAI-3D33. These models were supervised on population allele frequency data, a highly informative signal for pathogenicity prediction, which is also used to define labels in databases such as ClinVar34. Interestingly, although both AlphaMissense and PrimateAI-3D achieved top performance on ProteinGym and ClinVar, a simple post hoc adjustment of GPN-Star (V) predictions using gnomAD allele frequencies boosted its performance on ClinVar to surpass PrimateAI-3D and approach that of AlphaMissense (Fig. 2d and Supplementary Fig. 1).Pathogenic non-coding variantsWe next considered tasks involving non-coding variants, which are known to be especially challenging for predictive modelling. Here we demonstrate that GPN-Star is a powerful tool for identifying pathogenic non-coding variants in the human genome.In addition to previous genome-wide variant effect prediction (VEP) models, we included three prominent sequence-to-function models—Enformer35, Borzoi36 and AlphaGenome37—which are trained on extensive functional genomics data and have been widely applied to non-coding variant interpretation. We evaluated the models on classifying non-coding pathogenic variants from the OMIM38 and HGMD39 databases, both of which contain expert-curated annotations of human disease-associated variants. GPN-Star (M) achieved the best performance on both benchmarks (Fig. 2e,f). Notably, the sequence-to-function models performed substantially worse than alignment-based gLMs and phylogenetic models on this task, consistent with recent findings2,37. Stratified analysis indicated that GPN-Star consistently achieves top performance across variant types, with the largest advantage over sequence-to-function models in distal enhancers (Supplementary Figs. 2, 3).Given the critical role of promoter regions in transcription initiation and gene regulation, many specialized models have been developed to predict the impact of promoter variants. We therefore evaluated GPN-Star on promoter variants in OMIM and compared its performance with other methods, including three promoter-specific models: PromoterAI40, SpeciesLM41 and GPN-Promoter2. As shown in Fig. 2h, GPN-Star (M) demonstrated superior predictive performance compared with all competing models. In particular, the margin of improvement over the promoter-specific models is substantial.Fine-mapped variants in GWASGenome-wide association studies (GWAS) have been instrumental in identifying variants that contribute to genetic disease susceptibility. To further assess the use of GPN-Star, we evaluated its performance in classifying putatively causal versus non-causal missense variants resulting from fine-mapping of GWAS variants across 65 traits from the UK Biobank42. Among all competing models, GPN-Star (M) achieved the highest predictive performance on these fine-mapped missense variants. Notably, despite leveraging population allele frequency information, both AlphaMissense and PrimateAI-3D were substantially outperformed by GPN-Star in this task (Fig. 2c).We next evaluated the models on fine-mapped non-coding GWAS variants across 83 traits from the UK Biobank42, again assessing their ability to distinguish putatively causal versus non-causal variants. GPN-Star (M) continued to outperform all other models on this benchmark (Fig. 2g). Although sequence-to-function models (Enformer, Borzoi and AlphaGenome) exhibited moderate performance, Evo 2 showed relatively limited predictive value, as was previously observed2. Stratified analysis showed robust performance of GPN-Star across variant types (Supplementary Fig. 4). For fine-mapped variants located in promoter regions, GPN-Star (M) again outperformed all models, including the promoter-specific models PromoterAI, SpeciesLM and GPN-Promoter (Fig. 2h).Rare variant association testingHaving observed the strong performance of GPN-Star in predicting pathogenic and fine-mapped variants, we explored its use in rare variant association testing (RVAT), an important but challenging task in statistical genetics. The great interest in RVAT stems from the fact that most large-effect variants tend to be rare, as they are subject to negative selection pressure43. To overcome statistical power issues with low-frequency variants, RVAT is typically carried out at the gene level by aggregating variants in each gene44. These testing procedures often rely on variant annotations that reflect functional importance to prioritize variants in the aggregation. DeepRVAT5 is a recent method that uses a deep set network to integrate such variant annotations for RVAT. By combining an expressive deep learning framework with a powerful set of variant annotations, it demonstrated improved statistical power and computational efficiency compared with previous methods in extensive evaluations on UK Biobank whole-exome sequencing (WES) data5.We conducted RVAT experiments using DeepRVAT, enhancing the annotations used in the published version with all three GPN-Star predictions during both the training and association testing phases. Following the published benchmarking procedure5, we performed RVAT experiments using DeepRVAT on 34 quantitative traits in UK Biobank with the WES data from 161,822 unrelated individuals of European ancestry, retaining variants with minor allele frequency (MAF) <0.1% (Methods). To account for differences in stochastic model initialization, this procedure was run three times with different random seeds. As shown in Fig. 2j, adding the predictions from the GPN-Star models as variant annotations into DeepRVAT resulted in an increase in the number of discovered genes at a family-wise error rate <0.05 (on average across runs, 402 genes compared with 383 from the original method). Among the significant discoveries, there were also more gene–phenotype associations that replicated in conventional RVAT studies on UK Biobank with larger sample sizes45,46 (Methods) (on average 353 compared with 338), indicating high robustness of the further discoveries. Notably, the original DeepRVAT already includes state-of-the-art missense variant effect predictors such as AlphaMissense and PrimateAI, as well as the sequence-to-function model DeepSEA47, in the annotations. Nonetheless, GPN-Star seemed to offer complementary information and boosted the power of the tests.Complex trait heritabilityTo investigate the prediction of causal variants for complex trait GWAS, we turned to stratified linkage disequilibrium score regression (S-LDSC)48. S-LDSC is a principled way to estimate the informativeness of an annotation for complex trait heritability while leveraging the signal from all single nucleotide polymorphisms (SNPs), including those not confidently fine-mapped or genome-wide significant. It also serves as a backbone for functionally informed fine-mapping49 and polygenic risk scores50. Because S-LDSC requires scoring about 10 million variants, it was feasible to evaluate only the most scalable models (or those with precomputed scores). We binarized model scores to select a specific fraction of common variants. In our main analysis, we used the top 0.1% most constrained variants. We ran S-LDSC separately for each model score while conditioning on 96 baseline features and meta-analysed the results across 106 independent traits23.Since the early work on S-LDSC, conservation scores have been found to be the most enriched annotation for complex trait heritability48. More recently, primate-specific conservation—particularly PhastCons (P)—has emerged as the state-of-the-art for complex trait heritability enrichment7,8. Notably, GPN-Star (P) substantially improves on this result, followed by GPN-Star (M) (Fig. 3a). These improvements are even more striking when examining the heritability coefficient \({\tau }^{\star }\) (Fig. 3a), which quantifies the unique contribution of an annotation to heritability after adjusting for baseline features (Methods). Although per-trait enrichment estimates are noisier, GPN-Star (P) or GPN-Star (M) consistently rank at the top (Supplementary Fig. 5). This advance is particularly notable given the long-standing role of complex trait heritability enrichment as a meaningful benchmark in human genetics.Fig. 3: SNP heritability analyses in human complex traits.a, Informativeness of different models for complex trait heritability. We binarized model scores to select the top 0.1% most constrained common variants. We ran S-LDSC separately for each model annotation while conditioning on 96 baseline features and meta-analysed the results across 106 independent traits. Heritability enrichment is the proportion of heritability explained by an annotation divided by the size of the annotation. The conditional effect (\({\tau }^{\star }\)) measures the unique contribution to heritability that is not explained by existing annotations. b, Performance restricting model annotations to coding regions. c, Performance restricting model annotations to non-coding regions. d, Comparison varying the fraction of top constrained SNPs. e, Difference in enrichment between GPN-Star (P) and (M) as a function of estimated effective polygenicity across 27 traits. The P value is from a one-sided t-test on the Pearson correlation (d.f. = 25). Dashed line is an ordinary least squares fit. f, Types of common variants prioritized by GPN-Star (P). Variants are annotated by a combination of Ensembl consequences and ENCODE SCREEN candidate cis-regulatory elements. Only types with a proportion above 1% are shown. Odds ratio is given with respect to the bottom 99.9% of common variants. g, Comparison of tissue-agnostic and tissue-specific adaptations of GPN-Star (P), Enformer, Borzoi and tissue-specific baselines across selected trait groups (brain, 30 traits; blood/immune, 15 traits; liver, 15 traits). h, Comparison of tissue-agnostic and tissue-specific adaptations of GPN-Star (P), Enformer, Borzoi and tissue-specific baselines across selected individual traits. In a–d,g and h, annotations are defined on 9,997,231 reference SNPs, of which 5,961,159 are common. The bar and line heights are point estimates from S-LDSC. The error bars represent standard errors estimated from block jackknife over 200 genomic windows.The improvement of GPN-Star persists when restricting to top-ranked variants in coding and non-coding regions separately (Fig. 3b,c), with stronger gains in non-coding regions. We further performed an analysis of the best-performing model GPN-Star (P), stratifying the non-coding region by functional category, and found particularly strong improvements in enhancer regions (Extended Data Fig. 2). Moreover, across a range of binarization thresholds and all three evolutionary timescales, GPN-Star consistently outperforms previous conservation scores (Fig. 3d).Relevant evolutionary timescalesThe VEP results summarized in Fig. 2i underscore the importance of the timescale represented in the training data of evolutionary models. We further investigated how evolutionary timescale shapes the variants prioritized by GPN-Star models for complex trait heritability. We found that GPN-Star (P) preferentially captures variants contributing to highly polygenic traits51 (Extended Data Fig. 3), whereas GPN-Star (M) is relatively more informative for less polygenic traits, showing a connection between evolutionary timescale and the genetic architecture of complex traits (Fig. 3e). We also characterized the functional classes of variants prioritized by the top-performing model GPN-Star (P), finding large proportions of missense variants (33%, 170-fold enriched), followed by variants in regions with distal enhancer elements (26%, 2.4-fold enriched) and their flanks (8%, 0.35-fold depleted) (Fig. 3f). The importance of distal enhancer variants for complex trait heritability should not be understated, as they remain a major weakness of present sequence-to-function models52 (including the latest AlphaGenome37), as well as alignment-free gLMs (including the largest Evo 2; ref. 2). More analyses of model overlap and variant consequence enrichments are presented in Supplementary Information.Tissue-specific heritability signalsAs the final analysis of complex trait heritability, we investigated the performance of tissue-specific scores. Sequence-to-function models such as Enformer35 and Borzoi36 do not provide a single variant effect score but instead predict changes in activity across thousands of functional genomics tracks representing different assays and tissues. In our comparisons, we used Enformer and Borzoi scores aggregated across all tracks (tissue agnostic) as well as in seven specific tissues, following procedures from ref. 53 Given the low enrichment of tissue-agnostic Enformer and Borzoi scores when meta-analysed across all traits (Fig. 3a), we investigated the performance of tissue-specific scores. We meta-analysed these scores only across traits where we expected the tissue to be relevant (Supplementary Table 2). Tissue-specific Enformer and Borzoi scores consistently outperform their tissue-agnostic scores (Fig. 3g, Extended Data Fig. 4 and Supplementary Fig. 6).On the other hand, GPN-Star trained solely on DNA sequences is intrinsically tissue agnostic. Motivated by this observation, we devised a simple approach to incorporate tissue specificity into GPN-Star by only considering top-scoring variants located near tissue-specifically expressed genes (SEG) (inspired by LDSC-SEG54) or in tissue-specific candidate cis-regulatory elements (cCREs) (obtained from ENCODE v.453). Although this improved performance for most tissues—such as brain and blood/immune—it did not consistently surpass the performance of tissue-agnostic GPN-Star annotations, potentially because the present approach to incorporating tissue specificity remains suboptimal (Fig. 3g and Extended Data Fig. 4). Nevertheless, GPN-Star outperformed Enformer and Borzoi across all seven tissues. We also verified that the annotation by SEG and cCRE alone without model scores exhibited much lower heritability enrichment, confirming that the observed performance is still largely attributed to GPN-Star (Fig. 3g,h and Supplementary Fig. 6). Zooming in on example traits, brain-specific GPN-Star showed the highest enrichment for schizophrenia, blood/immune-specific GPN-Star performed best for lupus, and both tissue-agnostic and liver-specific GPN-Star yielded similar results for IGF1 (Fig. 3h). These findings highlight the importance of incorporating tissue-specific information in variant effect prediction.Interpretability of GPN-StarWe examined the representations learned by GPN-Star and found that its embeddings distinguish major genomic functional elements, with stronger separation for evolutionarily conserved regions (Supplementary Information; Supplementary Figs. 7–9). These results suggest that GPN-Star is aware of core functional elements of the genome when making predictions.Site-independent models such as PhyloP cannot, by definition, capture dependencies among nucleotides. By contrast, alignment-free gLMs have been shown to learn dependencies among well-known interacting elements, for example, in a transcription-factor binding site (TFBS) motif or between splice donors and acceptors41. To further probe understanding in GPN-Star of genomic syntax, we analysed learned nucleotide dependencies by systematically mutating each position and quantifying the resulting probability changes at other positions in the sequence41,55. This analysis showed biologically meaningful co-evolutionary signals, including dependencies among TFBSs, coding regions and splice sites, as well as evidence for primate-specific regulatory constraint. Representative examples are shown in Extended Data Fig. 5, with more case studies presented in Supplementary Figs. 10–13.These results demonstrate that GPN-Star can leverage co-evolutionary signals embedded in genomic sequence context to learn meaningful nucleotide dependencies that align with known functional dependencies, representing a notable advance over traditional conservation scores. Finally, we used nucleotide dependency to examine a few cases in our disease variant benchmarks where GPN-Star correctly predicts them to be deleterious, whereas PhyloP and PhastCons predict neutral scores (Extended Data Figs. 6, 7 and 8). We found that these variants do not exhibit strong site-wise conservation but reside in functionally important elements with distinct sequence contexts shown by nucleotide dependencies. As schematically depicted in Extended Data Fig. 9, GPN-Star effectively leverages these contextual signals to outperform traditional conservation metrics. For example, in the most extreme case of perfectly conserved sites, traditional conservation scores are unable to distinguish different variant types, whereas GPN-Star tends to predict more deleterious scores for higher impact variants (Supplementary Fig. 14).Genome-wide evolutionary constraintsTo more directly assess the connection between model predictions and evolutionary constraints in the genome, we leveraged allele frequency data from gnomAD v.3.1.2, which aggregates whole-genome sequencing samples from 76,156 human individuals30. Allele frequencies in the human population serve as informative indicators of selective constraint: more deleterious alleles tend to have lower frequencies as a result of purifying selection9,22,30.In this evaluation, we focused on comparison with PhyloP and PhastCons, which also learn evolutionary constraints from WGA data. To assess how well each model captures the relationship between allele frequency and constraint, we considered the vertebrate, mammal and primate versions of each model and obtained predictions for all gnomAD v.3 variants on chromosome 22, the held-out chromosome not used in training GPN-Star. We then compared the mean minor allele frequencies in several quantile bins defined by different models. As shown in Fig. 4a, across all three evolutionary timescales, variants in lower GPN-Star quantile bins have consistently lower average allele frequencies compared with those in corresponding PhyloP and PhastCons bins, suggesting that GPN-Star more accurately captures selective constraints in the human genome.Fig. 4: GPN-Star scores reflect evolutionary constraints on the human genome.a, Mean MAF for quantile bins ([0, 10−4], (10−4, 10−3], …, (10−1, 1]) defined by GPN-Star, PhyloP and PhastCons on the basis of the vertebrate, mammal and primate alignments at the gnomAD biallelic sites on chromosome 22. b,c, Enrichment of rare (singletons) versus common (MAF > 5%) gnomAD variants in the tail of deleterious scores (the threshold was chosen such that each score yielded 30 false discoveries): genome-wide enrichment (b) and enrichments stratified by molecular consequences of the variants (c). The rare variants were downsampled to match the number of common variants in each category. d, Performance comparison of GPN-Star with PhyloP and PhastCons at the three evolutionary timescales on correlation with Roulette62 mutation rate estimates on chromosome 22 (x axis) and correlation with Gnocchi30 constraint estimates (y axis). K, thousand; M, million.We next performed a more quantitative evaluation focusing on the most deleterious tails of the model score distributions, a regime especially important in many human genetics applications. We quantified the enrichment of rare variants relative to common variants in the most constrained tail predicted by each model. Here we defined rare variants as singletons and common variants as those with allele frequency >5%. Because rare variants are, on average, more deleterious, a model with more accurate constraint predictions should show higher enrichment of singletons in its most constrained tail. As shown in Fig. 4b and Supplementary Fig. 15a, all three GPN-Star models yielded substantially higher enrichment of rare variants than did PhyloP, PhastCons or CADD. Among the GPN-Star models, the vertebrate model showed the strongest enrichment overall and outperformed GPN-MSA (also trained on vertebrate genomes). When stratifying variants by molecular consequences, GPN-Star again achieved the highest enrichment in every category (Fig. 4c). Notably, GPN-Star (V) performed the best for missense variants, whereas GPN-Star (M) led for the synonymous and non-coding categories, mirroring trends observed in previous benchmarks. Lastly, our procedure for accounting for local mutation rate variation substantially reduced correlations with mutation rate estimates while improving performance across most downstream benchmarks (Fig. 4d and Extended Data Fig. 10; Supplementary Information).GPN-Star for model organismsGPN-Star requires only a multispecies WGA and the corresponding phylogeny to predict variant effects, making it readily applicable to species for which genome-wide variant interpretation has remained challenging. To demonstrate this, we trained GPN-Star models on five important model organisms: Mus musculus, Gallus gallus, Drosophila melanogaster, Caenorhabditis elegans and Arabidopsis thaliana. In each case, we showed that the learned constraints are highly informative for interpreting genetic variation in that species. For each of the five species, we collected a WGA dataset comprising aligned genomes from 18 to 135 species and applied training procedures similar to those used for the human models, with only minor adjustments (Supplementary Table 1).Owing to the scarcity of large-scale curated variant datasets with functional annotations in non-human species, we used population genetic data for evaluation. Specifically, we gathered population-level variant data for each of the five species (Methods) and examined the enrichment of rare variants in the most deleterious tail of the prediction distributions. Across all five species, GPN-Star scores exhibited substantially higher rare variant enrichment compared with PhyloP and PhastCons (Fig. 5a), indicating more accurate predictions of genome-wide evolutionary constraints. Across different variant categories based on molecular consequences, GPN-Star consistently showed the highest enrichment in almost all cases (Supplementary Figs. 15b–f and 16). In A. thaliana, the performance gap between GPN-Star and other models was smaller compared with the results for the other species, which could be due to either the small alignment or the known lower quality of WGAs in plants56.Fig. 5: Application of GPN-Star to non-human species.a, Evaluation of GPN-Star models for five non-human species on enrichment of rare versus common variants (thresholds in Methods) in five corresponding population genome databases in the tail of deleterious scores (the threshold was chosen such that each score yielded 30 false discoveries), compared with PhyloP and PhastCons fitted to the same alignments. b, Classification of MMRdb pathogenic variants versus WMGP common variants in M. musculus. c, Classification of FlyBase lethal variants versus DEST common variants in D. melanogaster. d, Classification of C. elegans lethal variants versus CaeNDR common variants. The performance metric used in b–d is the AUPRC. The bar height is the metric computed on the complete benchmark set and the error bars represent the 95% confidence intervals from 1,000 class-stratified bootstrap resamples. e, Nucleotide dependency map of the D. melanogaster model at a locus in the MSE enhancer annotated with known TFBS.For three of the species, we were able to gather curated pathogenic variants to carry out further evaluation. For the M. musculus model, we used pathogenic variants from the MMRdb database (Methods). GPN-Star outperformed both PhyloP and PhastCons in distinguishing these pathogenic variants from common variants in the population, both genome-wide (Fig. 5b) and in each variant category (Supplementary Fig. 16). Similarly, we collected experimentally validated lethal variants for D. melanogaster (from FlyBase57) and for C. elegans (from an experimental study58) and found that GPN-Star consistently outperformed the other models on both datasets (Fig. 5c,d and Supplementary Fig. 16). We also evaluated the mouse output head of Enformer35 on the non-coding variants in the mouse benchmarks (Supplementary Fig. 17). It was markedly outperformed by GPN-Star, echoing the non-coding VEP results for humans.Beyond variant prioritization, GPN-Star also enables exploration of learned functional elements and their co-evolution through nucleotide dependency analysis. As an example, we highlight learned nucleotide dependencies between TFBS at the well-characterized MSE enhancer locus in D. melanogaster (Fig. 5e). The coordinated activity of these transcription factors drives precise spatial patterning across the fly embryo59. This showcases how GPN-Star can serve as an efficient and low-cost tool for investigating functional elements and their dependencies across genomes, complementing and guiding experimental studies.DiscussionGPN-Star is a phylogeny-informed genomic language modelling framework that consistently outperforms existing methods for predicting functional constraint and deleterious variants across the human genome. Notably, GPN-Star achieved particularly strong predictive performance for distal enhancer variants, suggesting that evolutionary constraint information derived from cross-species comparisons can provide valuable signals for modelling this historically challenging class of regulatory element60. Beyond predictive accuracy, it learns biologically meaningful representations spanning functional elements from enhancers to TFBS and their co-evolutionary dependencies—all without any supervision. Classical WGA-based models, such as PhastCons18 and PhyloP19, have been indispensable tools to biologists and clinicians for over two decades. Our results indicate that GPN-Star complements and extends these tools, offering improved performance across diverse species and alignments.A central insight from our study is that the evolutionary timescale represented in the training data strongly influences the constraint learned by gLMs (Figs. 1c and 2i), corroborating and extending previous observations for classical phylogenetic methods7. Models trained on deeper timescales better capture coding and other highly constrained regions, whereas shallower timescales better inform rapidly evolving regulatory elements, including those constrained on primate-specific regulation8. The evolutionary timescale should therefore be considered an important design consideration for future gLMs.Like classical conservation scores, GPN-Star predictions quantify evolutionary constraint at specific timescales rather than general variant pathogenicity. Appropriate timescales should therefore be chosen according to the biological question. For integrating or selecting among models in predictive or discovery workflows, we recommend data-driven approaches that learn appropriate weights and transformations from task-specific data, as we demonstrated with DeepRVAT in our RVAT analysis.GPN-Star achieves state-of-the-art performance in functional constraint prediction with substantially smaller model and context sizes compared with existing single-sequence gLMs. This effectiveness probably stems from the explicit homology information provided by the alignment, which is not readily available to single-sequence models. Although we observed modest performance gains from increased model sizes, improvements from expanding context sizes were relatively small (Supplementary Fig. 18). A plausible explanation is that WGAs are constructed relative to a reference species and often consist of small, highly fragmented synteny blocks. As a result, when context size increases, nucleotides in a window may not be contiguous in the actual genomes of non-reference species, introducing misleading context information to the model. Exploring how to leverage larger context sizes with WGA data is a promising direction for future research. Advancement in model architecture design or alignment data structure is probably necessary.A practical limitation is that GPN-Star requires WGAs during inference. To facilitate its use, we have released genome-wide predictions through public repositories. Also, the reliance on a fixed alignment format makes the model less suitable for analysing indels and structural variants. On the other hand, single-sequence gLMs offer greater flexibility and applicability to a broad range of tasks and can naturally leverage large context and process indels and structural variants. We envision that future work can bridge these paradigms potentially through flexible homology retrieval mechanisms, allowing models to dynamically incorporate evolutionary information without relying on rigid alignment structures.Our results on pathogenicity prediction and complex trait heritability constitute a systematic evaluation of GPN-Star in the context of present knowledge of human genetic variation. However, the true value of the model depends on its capacity to enable biological discovery. We expect GPN-Star to hold great potential in human genetics applications by advancing our understanding of causal variants and human traits. Initial experiments demonstrate improved RVAT, and we expect even greater gains for non-coding analyses, in which variant prioritization remains challenging. GPN-Star may also improve functionally informed fine-mapping, polygenic risk prediction and phenotype prediction.Beyond identifying evolutionary constraints and deleterious variants, an important future direction is to study human-specific adaptations. The present GPN-Star models focus on capturing evolutionary signals across species, which lack the resolution for studying selection in human populations. We anticipate this to be an important direction for future work and would require learning from even more localized evolutionary data, such as human population genomes and archaic human genomes.It is remarkable how much can be learned from unlabelled DNA sequences alone. Evolutionary and functional genomics provide complementary views of genetic variation. Although evolutionary modelling outperformed present sequence-to-function approaches in our pathogenicity and heritability analyses, functional genomics remains essential for studying molecular phenotypes and tissue-specific regulation. Integrating these complementary sources of information into multimodal gLMs represents a promising direction.Large-scale efforts to sequence and align genomes across the tree of life are accelerating21,61. GPN-Star is well positioned to take advantage of this growth, as it scales readily to new alignments and may gain further capacity from increasingly diverse training data. In this age of rapidly expanding genomic data, we expect GPN-Star to be a valuable tool for leveraging these data to advance our understanding of genetic variation.MethodsSupplementary Information provides details of training data, model architecture, training setup, inference, mutation rate calibration, embedding analysis, nucleotide dependency analysis, evaluation benchmarks, heritability enrichment analysis and RVAT.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.