Three Major Key Results:

• Population Size & constraint
  Demonstrated that clades with large effective population sizes (Ne) show stronger constraint (lower dN/dS).
• Dynamic Selection Hotspots
  Proved that positive selection targets shift across clades (supporting epistasis).
• Taxon Power Multiplier
  Identified how dense species alignments act as power multipliers to recover previously missed selection.

Denser taxonomic databases introduce assembly and alignment errors. We classified errors into four categories:

• Kept Alignments (18,274 genes): Clean alignments passing all quality filters.
• Filtered Outliers (1,519 genes): Sequenced genomes containing extreme gaps or local misalignments.
• Filtered Pseudogenes (untracked in Kept): Non-functional regions with frameshift/inactivating errors that mimic selection under neutral drift.
• Too Few Sequences (160 genes): Less than 3 sequences, preventing phylogenetic models.

We set up an evolutionary association screen to detect genes where positive selection is enriched on terminal branches of echolocating lineages versus control lineages.

Methodological Details:

• Synonymous Control: Filtered out synonymous mutations (only non-synonymous selection tracked).
• Taxonomy Mapping: Integrated genus-based checks to prevent taxonomic misclassifications.
• Branch Control: Restricted selection calls to leaf branches (representing active, extant phenotypes).
• Statistical Test: Fisher's Exact Test per gene, with Benjamini-Hochberg FDR correction.

• Echolocators (110 spp): 411,386 positive selection events
• Background (632 spp): 1,858,010 positive selection events
• Echolocator-to-Background Ratio: ~1 to 4.52
• Note: The ratio reflects the larger branch subset of the background set, which is statistically normalized in our Fisher's Exact test background computations.

Running the screen against all 632 non-echolocating mammalian genomes controls for generic background rates and exposes specific molecular adaptations in echolocators:

Physiological Stress & Epigenetics:

• CD58 (Rank 1, Odds Ratio = 2.66): Key immune adhesion molecule. Selection suggests immune system remodeling, potentially to manage viral load in bats or diving-induced inflammatory stress in cetaceans.
• FAM111B (Rank 2, Odds Ratio = 2.46): Involved in DNA replication stress response. Echolocation creates high metabolic rates (flight/diving), producing systemic oxidative stress that requires enhanced DNA replication protection.
• NLRP13 (Rank 3, Odds Ratio = 2.12): Regulator of inflammatory responses, adapting echolocator tissues to survive hypoxia and transient tissue damage during intense activity.
• RSBN1 (Rank 5, Odds Ratio = 1.90): Epigenetic histone demethylase, indicating global transcriptional regulation changes.

Prestin (SLC26A5) is the canonical hearing gene mediating outer hair cell electromotility. Expanding sequence sampling to 742 mammal assemblies exposes key selective sites:

Prestin Selection Profile (Top 5 Sites):

• Codon 699: LRT = 22.06, p-value = 0.000007 (Significant)
• Codon 310: LRT = 17.04, p-value = 0.000084 (Significant)
• Codon 674: LRT = 15.46, p-value = 0.000186 (Significant)
• Codon 502: LRT = 13.47, p-value = 0.000507 (Significant)
• Codon 170: LRT = 13.28, p-value = 0.000560 (Significant)

We conducted a genome-wide Phylogenetic Generalized Least Squares (PGLS) screen to associate terminal branch selection events with simulated body mass (log grams) across species, controlling for phylogenetic covariance (shared ancestry).

Biological Interpretation of Top Hits:

• DAZAP1 (DAZ Associated Protein 1, Slope = 0.354): Crucial RNA-binding protein mediating spermatogenesis and growth dynamics. Positive selection on larger body size lineages suggests adaptation of development rate and reproductive lifespan scaling.
• LRRTM4 (Leucine Rich Repeat Transmembrane 4, Slope = 0.229): Synaptic cell adhesion molecule regulating neural wiring. Selection on larger body mass lines indicates structural neural adjustments during brain-body allometric scaling.
• GRM3 (Metabotropic Glutamate Receptor 3, Slope = -0.500): Selection associated with smaller body sizes, suggesting neuro-metabolic modifications to support higher relative mass-specific metabolic rates in smaller mammals.
• FGF13 (Fibroblast Growth Factor 13, Slope = 0.472): Direct role in microtubule stabilization and cardiac/neural development; adapts cellular growth processes to larger physical mass constraints.

A central result of this study is that scaling evolutionary analyses from sparse legacy alignments (about 40 species) to dense TOGA-level alignments (500 to 600 species) radically multiplies selection power.

Power Comparison Summary:

We identified high-confidence selection sites (classified as GOLD) that were completely lost in legacy datasets. Three representative case studies highlight this shift:

• ARID1B Site 498: Missing in 35-species tree vs. Highly Significant in 537-species tree.
• ARHGAP23 Site 3: Missing in 42-species tree vs. Highly Significant in 597-species tree.
• CABLES1 Site 554: Missing in 39-species tree vs. Highly Significant in 593-species tree.

The Invariant Site Mirage:

In the 35-species legacy alignment, this site exhibits exactly zero mutations. Because there is no variation, the evolutionary model has no data to distinguish purifying selection, a low mutation rate, or neutral drift. The likelihood surface is completely flat, leading to an LRT of 0.00.

By expanding the tree to 537 species in the TOGA dataset, the model uncovers the rich evolutionary history of this position, reconstructing 16 independent non-synonymous mutations and zero synonymous mutations. This high rate of purely amino-acid-altering change across the dense phylogeny provides the statistical replication necessary to confirm true selection with absolute confidence.

Distinguishing Drift from Selection:

In the legacy tree (42 species), the site has only a single non-synonymous mutation (occurring on the branch leading to the kangaroo rat Dipodomys ordii). A single mutation on a terminal branch is statistically indistinguishable from a random neutral substitution, yielding a low LRT of 0.30.

In the TOGA dataset (597 species), the site accumulates 17 independent non-synonymous mutations and zero synonymous mutations. Resolving these additional lineages provides the recurrent signals necessary to detect selection. The model estimates that 1.76% of branches undergo extremely intense positive selection (selection rate parameter beta_pos = 1858.37), resolving the signal with absolute confidence.

Widespread Selection Signal:

In the legacy alignment (39 species), this site has only a single non-synonymous mutation on an ancestral node. Similar to ARHGAP23, a single mutation provides insufficient evidence to rule out neutral evolution, leading to a p-value of 0.2902.

The dense TOGA tree (593 species) captures a massive accumulation of 31 independent non-synonymous mutations with zero synonymous changes. The model estimates that nearly the entire tree (99.0% of branches) evolves under diversifying selection (beta_pos = 3.95), indicating a widespread selection pressure that was completely lost under the sparse sampling of the legacy tree.

We evaluated which recently sequenced mammalian genomes (NCBI submission date on or after June 9, 2021) contribute the strongest selection signals to MEME results.

Top 5 Contributing Genomes:

1. Erinaceus europaeus (European hedgehog): Contributes 12,941 selection substitutions.
2. Ctenodactylus gundi (Gundi): Contributes 12,858 selection substitutions.
3. Typhlomys cinereus (Chinese pygmy dormouse): Contributes 11,191 selection substitutions.
4. Rhynchocyon petersi (Black and rufous elephant shrew): Contributes 9,894 selection substitutions.
5. Castor canadensis (American beaver): Contributes 9,204 selection substitutions.

Divergence & Bias:

• Primates/Cetacea branch lengths: Median is about 0.0009 substitutions/site.
• Rodentia/Chiroptera branch lengths: Median is about 0.006 substitutions/site.
• Unfiltered background analysis shows an artificial negative correlation (r is about -0.58) between terminal branch length and estimated background omega.

To resolve this math artifact, we implemented a branch-length filtering threshold, keeping only terminal branches with length greater than 0.003 substitutions/site.

After branch-length filtering (greater than 0.003), clades with larger effective population sizes (Ne) exhibit significantly stronger purifying selection (lower background omega):

Do evolutionary adaptive landscapes shift between deep mammalian lineages (epistasis), or remain constant across structural hotspots?

• Shared Selection Sites: Observed 91.84% is significantly lower than the expected 94.36% (Ratio = 0.973, p-value = 9.29e-201).
• Euarchontoglires-only Sites: Enriched by about 2-fold (Observed 4.05% vs Expected 2.14%).
• Laurasiatheria-only Sites: Enriched by 16% (Observed 3.97% vs Expected 3.42%).

Does adding more genomes to an alignment indefinitely yield linear returns in selection detection power?

• Immutable Housekeeping Genes: Genes annotated across all 700+ species represent core housekeeping machinery (e.g. ribosomes, translation, transcription).
• Tree Length (Mutational Depth): Evolutionary divergence (Tree Length) is the true predictor of selection detection power (r = 0.30, p-value = 3.05e-212), showing logarithmic return dynamics.

Taxon density challenges the paradigm of completely "invariant" negative controls. Housekeeping genes traditionally assumed to be under strict purifying selection across their entire sequence show site-specific episodic positive selection in our 742-species dataset:

• ACTB (Beta-actin): Standard structural component. Features 1 selection site under episodic positive selection (global dN/dS = 0.0033, P < 0.05).
• CALM1 (Calmodulin 1): Calcium sensor, 100% identity in standard genomes. Features 1 selection site (global dN/dS = 0.0022, P < 0.05).
• GAPDH: Standard expression reference. Features 3 selection sites (global dN/dS = 0.0402, P < 0.05).
• Ribosomal Subunits: RPL10, RPL19, RPL30, and RPS13 each contain 1 selection site under selection.

Truly Ultra-Conservative Genes (N = 476)

Only 4.6% of genes (476 / 10,247) remain completely devoid of selection (global dN/dS < 0.02, 0 selection sites, represented across ≥ 100 species). They are heavily enriched for:

• Cytoplasmic Translation (GO:0002181, Combined Score = 1484)
• mRNA Processing (GO:0006397, Combined Score = 168)
• Nucleosome Assembly (GO:0006334, Combined Score = 292)

Ultimate Negative Controls: Histone core proteins (e.g. H4C12, H3C10, H2BC4) remain 100% rigid across mammals, showing no codon variation.

Standard comparative genomic likelihood optimizations require hours per gene and fail to scale to VGP limits (thousands of species).

• Regularized Predictions: Rather than reproducing high-frequency asymptotic noise (such as alignment mismatches or tree length variations that inflate raw Likelihood Ratio Tests), the transformer acts as a biochemically regularized predictor.
• High Throughput: Bypasses raw numerical optimizations, accelerating selection screens by orders of magnitude to scale to VGP limits.