This is the results notebook for the dichromatism project — start here. It reports the central finding: sexual hair dichromatism has arisen and been lost repeatedly across the primate radiation, and where we can test its genetic basis, each origin is built through a different combination of two coupled gene modules — pigmentation and sex-hormone signaling.
It reads two things:
The phenotype — which primates are sexually dichromatic, scored across the radiation (§2), and where those gains and losses fall on the tree (§3).
Selection — per-gene episodic-selection scans on a curated 110-gene panel, run across the primate genomes on an HPC pipeline. The panel itself is justified in two upstream notebooks, which this one depends on but sits above in the reading order:
NB13 (13_sex_hormone_layer.ipynb) — why these sex-hormone genes, and how the hormone layer bridges into pigmentation.
NB14 (14_hormone_pigment_interface.ipynb) — the full 110-gene panel, its two modules, and the orthology screen behind every gene.
A reader who wants the finding needs only this notebook. A reader who wants why these genes follows the links down to NB13/NB14. The analysis flows the other way — the panel is built there and consumed here — but the science reads top-down from this page.
What this notebook establishes, in order:
§
Question
Result
2
How is dichromatism scored, and on how many species?
Coded across 238 species; 117 have sequenced genomes
Where selection can be tested, is one origin like another?
No — the powered origins use non-overlapping gene sets
5
Does an origin lean pigmentation or hormone?
Per-origin module balance, corrected for panel composition
6
Do independent origins reuse the same genes (convergence)?
No detectable shared signature — power-limited, but divergent where powered
0 — Input provenance
Every input below is a version-controlled file in this repository, traced to its producing notebook or pipeline and to its external source in DATA_SOURCES.md. The manifest is encoded as the INPUTS dict in the next cell so provenance is executable: the loader prints, for each file, what it is, where it came from, and a SHA256 of the exact bytes read. The phylogenetic results in §3 and §6 are recomputed here from the tree and coding — they are not read from any precomputed summary table.
primate time-tree (235 tips), pruned to the phenotype coding
grant-derived primate supertree; see DATA_SOURCES.md
config/primate_dichromatism_coding.csv
full-resolution dichromatism coding, 238 species
expert scoring of published descriptions/plates
analysis/data/dichromatism_coding.csv
the 117-species genome subset (dichromatic flag)
subset of the above with a sequenced genome
results/perorigin_v1/per_origin_K.csv
HPC per-origin RELAX (K) on the 110-gene panel
selection pipeline; panel justified in NB13/NB14
results/perorigin_v1/branch_rates.csv
HPC full-panel aBSREL per-branch selection
selection pipeline; panel justified in NB13/NB14
The phylo recompute is driven by nb15_phylo.R (committed beside this notebook); its outputs under dichromatism_synthesis/data/ and dichromatism_synthesis/figures/ are read below.
Show code
import os, hashlibimport pandas as pdREPO = os.environ.get("PIGNET_REPO", os.path.abspath(os.path.join(os.getcwd(), "..", "..", "..")))CG = os.path.join(REPO, "comparative-genomics")SYN = os.path.join(CG, "analysis", "dichromatism_synthesis") # this notebook's own data/figures# --- Provenance manifest: every input, where it comes from, how to regenerate it. ---INPUTS = {"tree": {"path": "comparative-genomics/analysis/coevolution_test/data/primate_phenotype_tree.nex","what": "Primate time-calibrated phylogeny (235 tips), the tree all rate/ASR estimates run on","source": "Published primate supertree pruned to the phenotype coding; see DATA_SOURCES.md","produced_by": "Upstream phylogenetics; this notebook only reads it", },"coding_full": {"path": "comparative-genomics/config/primate_dichromatism_coding.csv","what": "Full-resolution sexual-dichromatism coding, 238 species ""(hair_dichromatism_any + complete/state/type, natal coat, ontogenetic trajectory)","source": "Expert scoring of published species descriptions and plates","produced_by": "Phenotype coding effort; NB15 uses hair_dichromatism_any as the binary trait", },"coding_genome": {"path": "comparative-genomics/analysis/data/dichromatism_coding.csv","what": "The 117-species genome subset with the dichromatic flag ""(the species that also have a sequenced genome in the selection scan)","source": "Subset of coding_full restricted to species with a genome assembly","produced_by": "Genome-availability intersection; used to mark which origins are analysable", },"per_origin_K": {"path": "comparative-genomics/results/perorigin_v1/per_origin_K.csv","what": "Per-origin RELAX selection-intensity K per powered origin, on the pigmentation+""hormone panel (target 110 genes; the loaded row/gene count is printed at load ""time and grows as the clean-30 expansion and MYO5A/LYST giants land)","source": "HYPHY RELAX on codon alignments (HPC); panel justified in NB13/NB14","produced_by": "Selection pipeline (pulled from HPC); module+category columns from nb14_panel_justification.csv", },"branch_rates": {"path": "comparative-genomics/results/perorigin_v1/branch_rates.csv","what": "Full-panel aBSREL per-branch episodic-selection p-values on the 110-gene panel","source": "HYPHY aBSREL on codon alignments (HPC); panel justified in NB13/NB14","produced_by": "Selection pipeline (pulled from HPC)", },"relax_pooled": {"path": "comparative-genomics/results/perorigin_v1/relax_pooled_results.csv","what": "Pooled RELAX (all dichromatic origins as one foreground set) per gene, used as a ""QC cross-check on the per-origin fits: a per-origin hit with an extreme boundary ""K that is null when pooled is treated as a boundary artifact","source": "HYPHY RELAX, pooled foreground (HPC); panel justified in NB13/NB14","produced_by": "Selection pipeline (pulled from HPC)", },"sister_pairs": {"path": "comparative-genomics/analysis/module_selection/sister_pairs.csv","what": "Each dichromatic taxon paired with its closest monochromatic relative on the tree","source": "Nearest-monochromatic-sister search on the phenotype tree","produced_by": "analysis/module_selection (sister-pair construction)", },"gene_modules": {"path": "comparative-genomics/analysis/module_selection/data/gene_modules.csv","what": "Module label (pigmentation / hormone) for every gene in the 107-gene panel","source": "Panel definition (NB13/NB14)","produced_by": "analysis/module_selection (panel module map)", },# --- Pre-built figures (aBSREL tip-level module views over ALL origins, incl. single-species) ---# These four are rendered by scripts in comparative-genomics/analysis/module_selection/ from the# SAME final panel this notebook loads (branch_rates.csv, 105 genes; module column, NOT set).# They are read here as image inputs, exactly as NB15 reads NB13/NB14 figures."fig_module_balance": {"path": "comparative-genomics/analysis/dichromatism_synthesis/figures/nb15_fig_module_balance.png","what": "Per-origin module balance (aBSREL tip-level, all 14 origins) as a diverging bar chart","source": "aBSREL episodic selection on the 107-gene panel (HPC)","produced_by": "comparative-genomics/analysis/module_selection/figures.py", },"fig_circular_tree": {"path": "comparative-genomics/analysis/dichromatism_synthesis/figures/nb15_fig_circular_tree_balance.png","what": "The same module-balance metric painted on the circular primate tree (per-tip)","source": "aBSREL episodic selection on the 107-gene panel (HPC)","produced_by": "comparative-genomics/analysis/module_selection/tree_figures.R", },"fig_sister_contrast": {"path": "comparative-genomics/analysis/dichromatism_synthesis/figures/nb15_fig_sister_pair_contrast.png","what": "Selection in each dichromatic taxon vs its closest monochromatic sister, gene x taxon matrix","source": "aBSREL episodic selection (HPC); sister pairs from the phenotype tree","produced_by": "comparative-genomics/analysis/module_selection/sister_figures.py", },"fig_sister_network": {"path": "comparative-genomics/analysis/dichromatism_synthesis/figures/nb15_fig_sister_network_diff.png","what": "Where the dichromatic-vs-sister selection difference sits in the coupled gene network","source": "aBSREL episodic selection (HPC); STRING v12 human edges (score>=0.4)","produced_by": "comparative-genomics/analysis/module_selection/sister_figures.py", },}def _sha(path):return hashlib.sha256(open(path, "rb").read()).hexdigest()[:16]def load_input(key, reader=pd.read_csv):"""Load a manifest input, printing provenance + a checksum of the exact bytes read. Returns None (with a clear notice) if the file is not present yet — the HPC selection tables land after the tree-only sections are built.""" m = INPUTS[key]; fp = os.path.join(REPO, m["path"])ifnot os.path.exists(fp):print(f"[{key}] {m['path']}\n NOT PRESENT YET — {m['what']}\n"f" (produced by: {m['produced_by']})\n")returnNoneprint(f"[{key}] {m['path']}")print(f" what : {m['what']}")print(f" source : {m['source']}")print(f" produced : {m['produced_by']}")if key =="tree":print(f" loaded : Nexus tree | sha256[:16]={_sha(fp)}\n");return fp df = reader(fp)print(f" loaded : {df.shape[0]} rows x {df.shape[1]} cols | sha256[:16]={_sha(fp)}\n")return dfTREE_PATH = load_input("tree")CODING = load_input("coding_full")GENOME = load_input("coding_genome")
[tree] comparative-genomics/analysis/coevolution_test/data/primate_phenotype_tree.nex
what : Primate time-calibrated phylogeny (235 tips), the tree all rate/ASR estimates run on
source : Published primate supertree pruned to the phenotype coding; see DATA_SOURCES.md
produced : Upstream phylogenetics; this notebook only reads it
loaded : Nexus tree | sha256[:16]=cfed56b5699eb900
[coding_full] comparative-genomics/config/primate_dichromatism_coding.csv
what : Full-resolution sexual-dichromatism coding, 238 species (hair_dichromatism_any + complete/state/type, natal coat, ontogenetic trajectory)
source : Expert scoring of published species descriptions and plates
produced : Phenotype coding effort; NB15 uses hair_dichromatism_any as the binary trait
loaded : 238 rows x 12 cols | sha256[:16]=96c6d7e94e02bdae
[coding_genome] comparative-genomics/analysis/data/dichromatism_coding.csv
what : The 117-species genome subset with the dichromatic flag (the species that also have a sequenced genome in the selection scan)
source : Subset of coding_full restricted to species with a genome assembly
produced : Genome-availability intersection; used to mark which origins are analysable
loaded : 117 rows x 4 cols | sha256[:16]=9f064b9d9a6c7c3a
1 — The headline
Sexual dichromatism — males and females differing in hair or pelage color — is scattered across the primate tree, not clustered in one clade. Coded across 238 species, it appears in lemurs (Eulemur), colobines (Trachypithecus, Colobus), gibbons (Nomascus, Hylobates), and a handful of other lineages, separated by long monochromatic stretches. That distribution is the starting observation: whatever builds dichromatism has been assembled repeatedly and independently, and lost repeatedly too.
Two questions, two kinds of data. The phylogeny answers where, when, and how often the trait arose — that it is labile, lost far faster than gained, and assembled ~15 times (§3). None of that needs sequence data. The reason we sequenced is a question the tree cannot answer: what molecular change happened when dichromatism arose — and here the answer is genuinely surprising. In birds, sexual dichromatism repeatedly evolves through MC1R, a single large-effect coat-color gene. In primates, MC1R is not the recurrent route: across the panel it shows no pooled selection signal (K = 0.45, n.s.), surfacing in just one lineage (§4). There is no single “dichromatism gene.” Instead, the origins we can test recruited different gene sets from a coupled pigmentation–hormone system, with no shared molecular signature where the comparison is powered (§4–§6). That is the payoff of the genetics: not the pattern of origins (the tree gives that), but the finding that a trait built one way in birds is built many different ways in primates.
Roadmap: §3 quantifies the lability and origin count from the tree; §4–§5 test which genes and which module shifted at the origins with enough sequenced species (only two carry a signal); §6 asks whether independent origins converged on the same genes (they do not, where powered). Whether origins differ systematically in architecture is left as an open, underpowered question — only two of ~15–19 origins carry a per-origin signal — not asserted.
2 — Phenotype: coding and provenance
Two codings, one trait. Dichromatism is scored at full taxonomic resolution across 238 species (coding_full), using hair_dichromatism_any as the binary trait (1 = males and females differ in hair color). A 117-species subset (coding_genome) marks the species that also have a sequenced genome — these are the ones a selection scan can actually run on. The two codings agree on the phenotype; they differ only in taxonomic coverage, and that difference matters for one number (the origin count, §3), so both scopes are reported throughout rather than silently merged.
The tree (tree, 235 tips) and the coding are joined on the binomial with underscores substituted for spaces; the intersection used for all rate/ASR work is reported below.
Show code
# Join tree and coding; report the overlap used for all downstream phylo work.import redef _tips_from_nexus(path): txt =open(path).read()returnset(re.findall(r"[A-Z][a-z]+_[a-z]+", txt))if TREE_PATH isnotNoneand CODING isnotNone: tree_tips = _tips_from_nexus(TREE_PATH) CODING = CODING.copy() CODING["tip"] = CODING["species_binom"].str.replace(" ", "_", regex=False) coded = CODING.dropna(subset=["hair_dichromatism_any"]) overlap =sorted(set(tree_tips) &set(coded["tip"])) n_dich =int(coded.set_index("tip").loc[overlap, "hair_dichromatism_any"].astype(int).sum()) n238 =len(CODING); n117 =0if GENOME isNoneelselen(GENOME)print(f"full coding (coding_full) : {n238} species")print(f"genome subset (coding_genome) : {n117} species")print(f"tree tips : {len(tree_tips)}")print(f"tree n coding (analysis set) : {len(overlap)} tips, {n_dich} dichromatic")print("\nAll rate and ancestral-state estimates below run on the "f"{len(overlap)}-tip intersection.")
full coding (coding_full) : 238 species
genome subset (coding_genome) : 117 species
tree tips : 235
tree n coding (analysis set) : 224 tips, 25 dichromatic
All rate and ancestral-state estimates below run on the 224-tip intersection.
3 — The trait is evolutionarily labile
Method. A two-state Markov model (Mk) is fit to hair_dichromatism_any on the intersection tree, comparing equal rates (ER: one rate for gain and loss) against all-rates-different (ARD: separate gain and loss rates). Model fit is compared by AIC; the origin count is read from stochastic character maps under the preferred model, and cross-checked against a model-free count of maximal dichromatic clades. All of this is recomputed by nb15_phylo.R; its result tables are read below.
Table 1. Model comparison for how dichromatism changes over evolutionary time. ER = one rate for gains and losses; ARD = separate gain and loss rates. Lower AIC = better-fitting model; a gap (ΔAIC) above ~4 is decisive. ARD wins by ΔAIC ≈ 20, and its two rates show that the trait is lost far faster than it is gained (loss:gain ≈ 9×). The origin count — how many separate times dichromatism appeared — is read from simulations under the winning model and reported for both coding scopes, because that one number depends on how many species are included.
Show code
# Read the recomputed lability fits + origin estimates (from nb15_phylo.R).def _read_syn(name): fp = os.path.join(SYN, "data", name)return pd.read_csv(fp) if os.path.exists(fp) elseNoneFITS = _read_syn("nb15_lability_fits.csv")ORIGINS = _read_syn("nb15_origin_estimates.csv")if FITS isnotNone: ard = FITS[FITS.model =="ARD"].iloc[0]print("Mk model comparison (delta = ER - ARD):")for _, r in FITS.iterrows():print(f" {r['model']:3s} logLik={r['logLik']:8.2f} AIC={r['AIC']:8.2f}")print(f"\nARD preferred by dAIC = {ard['dAIC_ER_ARD']:.2f}")print(f" gain rate (0->1) = {ard['gain_0to1']:.4f}")print(f" loss rate (1->0) = {ard['loss_1to0']:.4f}")print(f" loss : gain = {ard['loss_gain']:.1f}x (losses far outpace gains)")if ORIGINS isnotNone:print("\nOrigin count - reported under both coding scopes (differ by sampling, not phenotype):")for _, r in ORIGINS.iterrows(): v =""if pd.isna(r["n_origins"]) elsef"{r['n_origins']:.0f}"print(f" {r['scope']:38s}{v:>4s} [{r['method']}]")
Mk model comparison (delta = ER - ARD):
ER logLik= -80.25 AIC= 162.49
ARD logLik= -69.35 AIC= 142.70
ARD preferred by dAIC = 19.79
gain rate (0->1) = 0.0273
loss rate (1->0) = 0.2484
loss : gain = 9.1x (losses far outpace gains)
Origin count - reported under both coding scopes (differ by sampling, not phenotype):
full coding (238 species) 19 [maximal dichromatic clades (topology, model-free)]
full coding, genome-sampled clades 18 [clades with >=1 genome-sampled tip]
genome-subset coding (117 species) 15 [maximal dichromatic clades (topology, model-free)]
stochastic map ARD (mean) 31 [500 simmaps (flicker-inflated by high loss rate)]
Result. The ARD model is strongly preferred (ΔAIC ≈ 19.8): a single-rate model is a poor description. The asymmetry is large and in one direction — dichromatism is lost roughly nine times faster than it is gained. A trait that is hard to gain but easy to lose, arising in scattered clades, is the signature of a labile, repeatedly-assembled phenotype rather than one ancient origin retained.
Origin count — an estimator- and sampling-dependent number. How many times dichromatism arose depends on both the coding scope and the method:
On the 117-species genome subset — the coding the selection scan actually runs on — a model-free count of maximal dichromatic clades gives ~15 origins. This is the figure used in the project’s README and walkthrough.
On the full 238-species coding, the finer taxonomic sampling resolves 19 distinct dichromatic clades (18 containing at least one genome-sampled species). The extra origins are real clades the genome subset simply did not sample.
Model-based counts bracket these: ML ancestral-state reconstruction under ARD collapses the estimate (it paints dichromatism deep and loses it repeatedly — a known artifact of the high loss rate, and the reason a naïve “4 deep origins” reading is wrong), while stochastic mapping inflates it (mean ≈ 31, from transient 0→1→0 flicker on short branches).
The defensible statement is therefore a range, not a point: on the order of 15–19 independent origins, ~14–18 of them with a sequenced genome. The exact value is sensitive to how densely the clade is sampled — which is itself part of the lability story.
Reading the density map without circularity
A fair objection to the map below (Fig 3): it is a reconstruction fitted to the observed tip states, so at the tips the posterior is pinned to what we coded — a red-painted branch ending in a red observed dot is close to tautological, and proves nothing on its own. The map carries real information in two places that are not circular: (i) at internal branches, where it infers where gains and losses happened; and (ii) in flagging species whose state is unexpected given their relatives.
We test (ii) directly with a leave-one-out check: predict each species’s state from the rest of the tree — its parent node’s ancestral reconstruction propagated down its own branch — and compare to what is observed. This asks each species to be explained by its neighbors, not by itself, so it cannot be circular. nb15_phylo.R computes it (nb15_loo_surprise.csv); the most surprising species are shown in Table 2.
One caveat stated plainly: because the fitted model already says gains are rare and recent (loss:gain ≈ 9×), it will call most gains “surprising” by construction — so the count of surprising gains is not independent of the lability result, it is the same signal per species. What the test adds is which species, and the asymmetry (many surprising gains, almost no surprising losses) as an independent readout of the same high-loss dynamic. A fully out-of-sample version — predicting primates not in our sample and scoring against their known coats — is the clean next step and is noted in the limitations.
Show code
# Load the leave-one-out surprise table (computed reproducibly in nb15_phylo.R), join the# developmental coding from the same CSV that defines the trait, and show the species whose# observed state most disagrees with the prediction from the rest of the tree._surf = os.path.join(SYN, "data", "nb15_loo_surprise.csv")if _surf and os.path.exists(_surf): SUR = pd.read_csv(_surf)# bring in natal-coat / ontogeny coding from the trait CSV (keyed on the same binomial->tip) _cod = CODING.copy() _dev = _cod.set_index("tip")[["natal_coat", "ontogenetic_trajectory_color"]] SUR = SUR.merge(_dev, left_on="tip", right_index=True, how="left") _lab = {1: "dichromatic", 0: "monochromatic"} top = SUR.head(12).copy() top["observed"] = top["observed"].map(_lab)print("Most phylogenetically surprising species, with developmental coding:\n")print(top[["tip", "observed", "predicted", "surprise","natal_coat", "ontogenetic_trajectory_color"]].to_string(index=False)) ng =int(((SUR.observed ==1) & (SUR.predicted <0.5)).sum()) nd =int((SUR.observed ==1).sum()) nl =int(((SUR.observed ==0) & (SUR.predicted >0.5)).sum())print(f"\nUnexpected dichromatic origins (dichromatic, predicted < 0.5): {ng} of {nd}")print(f"Reversal candidates (monochromatic, predicted > 0.5): {nl}")# near-total coupling of natal coat with adult dichromatism (largely definitional; see prose) dich = SUR[SUR.observed ==1]; mono = SUR[SUR.observed ==0]print(f"\nNatal coat present: {int((dich.natal_coat ==1).sum())}/{len(dich)} dichromatic vs "f"{int((mono.natal_coat ==1).sum())}/{len(mono)} monochromatic species.")print("Ontogenetic trajectory among the surprising dichromatic origins:")print(dich[dich.predicted <0.5]["ontogenetic_trajectory_color"] .value_counts().to_string())else:print("nb15_loo_surprise.csv not built yet - run nb15_phylo.R")
Table 2. The phylogenetically surprising species (leave-one-out), with developmental coding. From nb15_loo_surprise.csv joined to the trait CSV. observed = the coded state; predicted = P(dichromatic) estimated from the rest of the tree (parent-node reconstruction propagated down the species’ own branch); surprise = |observed − predicted|, so 1.0 = maximally unexpected. natal_coat = whether infants wear a distinct natal coat; ontogenetic_trajectory_color = how the coat matures (bidirectional_maturation = infant differs from both adult sexes and matures toward each; male_maturation = infant resembles the adult female, males diverge at maturity; early_dimorphism = sexes differ from early on). The top of the list is dominated by dichromatic species sitting in otherwise-monochromatic genera (Macaca arctoides, Alouatta caraya, Erythrocebus patas, Pithecia pithecia, …) — the single-species independent origins that make dichromatism look scattered. Only one clear reversal appears (a monochromatic species predicted dichromatic): Trachypithecus delacouri, embedded in the dichromatic langur radiation.
The surprising origins are whole developmental programs, not just adult color
The surprising dichromatic origins do not merely have differently-colored adult males and females — they carry a full ontogenetic program: a distinct infant natal coat that then matures into the sex-specific adult pattern. Almost every dichromatic species in the tree carries a natal coat (24 of 25) versus ~6% of monochromatic species (11 of 199), and the surprising and expected dichromatic species carry it at the same near-universal rate.
A definitional caveat, stated plainly. This coupling is partly built into the coding: if the adults are dichromatic and the species is not dichromatic from birth, an infant coat distinct from at least one adult sex must exist — so “has a natal coat” is close to a logical consequence of “adults are dichromatic,” not an independent correlate. What is not definitional, and is the point here, is the trajectory: the surprising origins are mostly bidirectional_maturation and male_maturation, i.e. dichromatism is reached by sex-specific maturation away from a shared juvenile coat, not by males and females differing from birth. That means each of these scattered origins switched on the same kind of developmental module — a coordinated natal-coat-to-adult trajectory — rather than an ad-hoc adult-color change. It reframes the lability finding: what is gained and lost ~15 times is a developmental program, not a paint job.
Figure 3 — Where gains and losses fall (stochastic density map)
The figure below paints, on the primate tree, how likely dichromatism was along each branch — so gains and losses of the trait are visible as where red appears and disappears. Full caption follows the figure.
Show code
from IPython.display import Image, display_dm = os.path.join(SYN, "figures", "nb15_densitymap.png")if os.path.exists(_dm): display(Image(filename=_dm))else:print("nb15_densitymap.png not built yet - run nb15_phylo.R")
Figure 3. Where on the primate tree dichromatism was gained and lost (density map).
What the data are. The tree is the 224-species primate phylogeny used throughout this notebook (the intersection of the 235-tip phylogeny and our 238-species dichromatism coding, §2); 25 of those species are coded dichromatic. The trait is hair_dichromatism_any (males and females differ in hair color: yes/no).
What a density map is. Every branch of the tree is colored by how likely it is that dichromatism was present along that branch. To estimate this we simulate the trait’s evolution down the tree 500 times (stochastic character maps under the best-fitting ARD model from §3) and average the results, so a branch that comes out dichromatic in most simulations is painted red and one that is rarely dichromatic stays blue.
Branch color runs from blue (probability ≈ 0, monochromatic) through magenta to red (probability ≈ 1, dichromatic); the scale bar is at lower-left. A red dot at a tip marks a living species that is actually observed / coded dichromatic — this is the reality check on the reconstruction: the painted-red branches should end in red-dotted tips, and they do, so the model is tracking the real states rather than inventing signal. The colored strip on the right labels the major primate groups for orientation (Old World monkeys, apes, gibbons, New World monkeys, tarsiers, lemurs, lorises & galagos). Individual species names are left off so the 224-species tree stays legible (§2).
Red concentrates in short, terminal patches — the Trachypithecus langurs within the Old World monkeys, the Nomascus/Hylobates gibbons, and Eulemur among the lemurs — separated by long blue (monochromatic) internodes. That scattering, rather than one deep red clade, is the visual signature of the many-independent-origins, high-loss pattern quantified in §3.
On circularity. At the tips this map is pinned to the observed states, so the tip-level match is expected by construction and is not evidence on its own (see “Reading the density map without circularity” above). The non-circular readout — which species are unexpected given their relatives — is in Table 2: the red patches are overwhelmingly single-species gains their neighbors would not predict, with essentially one reversal.
4 — Per-origin architecture
The question. Of the independent origins, only those with ≥2 sequenced dichromatic tips carry enough branches for a per-origin RELAX test — three do: Trachypithecus (origin 7), Nomascus (origin 8), and Eulemur (origin 14). For each, which panel genes show a shift in selection intensity (K) along the origin’s branches?
The result (the loaded panel size is printed below and grows toward 110 as the clean-30 pigmentation expansion lands — the two giant genes MYO5A/LYST arrive last): the three powered origins do not share a gene set. Trachypithecus shows a multi-gene, both-module signal (including a lineage-specific MC1R shift); Nomascus concentrates on a small pigmentation set — POMC, HGF, HRAS — of which POMC (§5b) sits at the pigmentation–hormone interface; Eulemur shows no gene passing the per-origin threshold. Same phenotype, different genetic routes — quantified in §6. The gene lists are printed below; the fig_per_lineage_genes panel is regenerated against the completed 110-gene tables.
Show code
PER_ORIGIN = load_input("per_origin_K")BRANCH = load_input("branch_rates")if PER_ORIGIN isnotNone: origins =sorted(PER_ORIGIN['origin_id'].unique()) if'origin_id'in PER_ORIGIN else'?'print(f"per_origin_K: {PER_ORIGIN.shape[0]} rows, "f"{PER_ORIGIN['gene'].nunique()} genes, origins {origins}")if"p_BH"in PER_ORIGIN.columns: sig = PER_ORIGIN[PER_ORIGIN.p_BH <0.05]for o in origins: g =sorted(sig[sig.origin_id == o].gene)print(f" {o}: {len(g)} genes p_BH<0.05 - {g if g else'(none)'}")else:print("Awaiting per-origin tables (cluster appending clean-30 into perorigin_v1).")
[per_origin_K] comparative-genomics/results/perorigin_v1/per_origin_K.csv
what : Per-origin RELAX selection-intensity K per powered origin, on the pigmentation+hormone panel (target 110 genes; the loaded row/gene count is printed at load time and grows as the clean-30 expansion and MYO5A/LYST giants land)
source : HYPHY RELAX on codon alignments (HPC); panel justified in NB13/NB14
produced : Selection pipeline (pulled from HPC); module+category columns from nb14_panel_justification.csv
loaded : 273 rows x 13 cols | sha256[:16]=03afa4def0c6306a
[branch_rates] comparative-genomics/results/perorigin_v1/branch_rates.csv
what : Full-panel aBSREL per-branch episodic-selection p-values on the 110-gene panel
source : HYPHY aBSREL on codon alignments (HPC); panel justified in NB13/NB14
produced : Selection pipeline (pulled from HPC)
loaded : 12806 rows x 13 cols | sha256[:16]=f1d228d89f0e3016
per_origin_K: 273 rows, 106 genes, origins ['origin_14', 'origin_7', 'origin_8']
origin_14: 0 genes p_BH<0.05 - (none)
origin_7: 13 genes p_BH<0.05 - ['ASIP', 'EDN3', 'HPS4', 'HSD17B7', 'KISS1R', 'MC1R', 'SCARB1', 'SHBG', 'SPTLC2', 'SRD5A1', 'SRD5A3', 'TFAP2A', 'TYR']
origin_8: 3 genes p_BH<0.05 - ['HGF', 'HRAS', 'POMC']
Figure 4 — Per-origin selection architecture
For each origin of dichromatism that contains enough sequenced species to test, the figure below shows which panel genes shifted in selection and in which module. Full caption follows the figure.
Show code
# QC gate (from the cluster's pooled-RELAX evaluation): a per-origin hit with an extreme K# (K>20 or K<0.05) that is NULL in the pooled analysis is a boundary artifact and is dropped# (e.g. HPS4: per-origin K=30 but pooled K=1.05, p_BH=1.0). Bars are plotted at the pooled K where# a pooled estimate exists (so a boundary K is shown at its de-inflated value), and marked by# whether the pooled analysis corroborates the hit.import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.patches import Patch_arch = os.path.join(SYN, "figures", "nb15_per_origin_arch.png")POOL = load_input("relax_pooled")if PER_ORIGIN isnotNoneand"p_BH"in PER_ORIGIN.columns: _labels = {"origin_7": "Trachypithecus\n(origin 7, 8 tips)","origin_8": "Nomascus\n(origin 8, 3 tips)","origin_14": "Eulemur\n(origin 14, 2 tips)"} _sig = PER_ORIGIN[PER_ORIGIN.p_BH <0.05].copy() _gmod = PER_ORIGIN.drop_duplicates("gene").set_index("gene")["module"].to_dict() _mc = {"pigmentation": "#c0662e", "hormone": "#3b6ea5"}# pooled lookup for QC + de-inflation _pl = POOL.set_index("gene")[["K", "p_BH"]].to_dict("index") if POOL isnotNoneelse {}def _keep(row): # drop boundary-K hits that pooled analysis calls null k = row["K"]; pk = _pl.get(row["gene"], {}) boundary = (k >20) or (k <0.05) pooled_null = pk.get("p_BH", 0) >0.05returnnot (boundary and pooled_null) _sig = _sig[_sig.apply(_keep, axis=1)].copy()# plot K: use pooled K when the per-origin K is boundary-inflated and a pooled value existsdef _plotK(row): k = row["K"]; pk = _pl.get(row["gene"], {})if (k >20or k <0.05) and"K"in pk and pk["K"] == pk["K"]:return pk["K"]return k _sig["plotK"] = _sig.apply(_plotK, axis=1) _sig["corrob"] = _sig["gene"].map( # pooled corroborates the intensification/relaxation?lambda g: _pl.get(g, {}).get("p_BH", 1.0) <0.05) _origs = [o for o in _labels if o in PER_ORIGIN.origin_id.values] fig, axes = plt.subplots(1, len(_origs), figsize=(4.5*len(_origs), 5.4), gridspec_kw={"wspace": 0.6})iflen(_origs) ==1: axes = [axes]for ax, o inzip(axes, _origs): s = _sig[_sig.origin_id == o].copy()iflen(s) ==0: ax.text(0.5, 0.5, "no gene passes\np(BH) < 0.05", ha="center", va="center", transform=ax.transAxes, fontsize=12, color="#7f8c8d", style="italic") ax.set_title(_labels[o], fontsize=12, fontweight="bold"); ax.axis("off");continue s["logK"] = np.log2(s["plotK"].clip(lower=0.05)) s = s.sort_values(["module", "logK"]) y = np.arange(len(s))# filled bar = pooled-corroborated; hatched/pale = origin-specific (not corroborated pooled)for yi, (_, r) inzip(y, s.iterrows()): col = _mc[_gmod.get(r.gene, "pigmentation")] ax.barh(yi, r["logK"], color=col if r["corrob"] else"white", edgecolor=col, linewidth=1.6, height=0.7, hatch=Noneif r["corrob"] else"///") ax.set_yticks(y); ax.set_yticklabels(s.gene, fontsize=9, fontstyle="italic") ax.axvline(0, color="black", lw=0.8) ax.set_title(_labels[o], fontsize=12, fontweight="bold") ax.set_xlabel("log$_2$ K (selection intensity)", fontsize=9)# directional cues placed on the side they describe: relaxed to the LEFT, intensified RIGHT ax.annotate("← relaxed", xy=(0, -0.14), xycoords="axes fraction", ha="left", va="top", fontsize=8, color="#555") ax.annotate("intensified →", xy=(1, -0.14), xycoords="axes fraction", ha="right", va="top", fontsize=8, color="#555") ax.spines[["top", "right"]].set_visible(False) fig.suptitle("Per-origin selection architecture: which genes shift, in which module", fontsize=13, fontweight="bold", y=1.03) fig.legend(handles=[Patch(facecolor=_mc["pigmentation"], edgecolor="black", label="pigmentation"), Patch(facecolor=_mc["hormone"], edgecolor="black", label="hormone"), Patch(facecolor="white", edgecolor="#555", hatch="///", label="not corroborated by pooled RELAX (origin-specific)")], loc="lower center", ncol=3, frameon=False, bbox_to_anchor=(0.5, -0.10), fontsize=9) fig.savefig(_arch, dpi=150, bbox_inches="tight") plt.show()print("dropped as boundary+pooled-null:",sorted(set(PER_ORIGIN[PER_ORIGIN.p_BH <0.05].gene) -set(_sig.gene)) or"none")else:print("Per-origin architecture figure builds when the tables are present.")
[relax_pooled] comparative-genomics/results/perorigin_v1/relax_pooled_results.csv
what : Pooled RELAX (all dichromatic origins as one foreground set) per gene, used as a QC cross-check on the per-origin fits: a per-origin hit with an extreme boundary K that is null when pooled is treated as a boundary artifact
source : HYPHY RELAX, pooled foreground (HPC); panel justified in NB13/NB14
produced : Selection pipeline (pulled from HPC)
loaded : 109 rows x 11 cols | sha256[:16]=48dca781e1be583f
dropped as boundary+pooled-null: ['HPS4']
Figure 4. Genes under selection at each independent origin of dichromatism.
What the data are. We scanned a curated panel of ~110 pigmentation and sex-hormone genes (built and justified in NB14) for signatures of natural selection, using coding sequences extracted from 117 published primate genomes. Selection was tested with RELAX (HyPhy), which asks whether selection on a gene became stronger (intensified) or weaker (relaxed) along a specified set of branches; its statistic K summarizes that shift, and we plot log₂ K so that intensification (K > 1) and relaxation (K < 1) are symmetric about zero. Values are read from results/perorigin_v1/per_origin_K.csv (per-origin fits) cross-checked against relax_pooled_results.csv (all origins pooled).
What the panels are. Dichromatism arose ~15 times (§3), but only origins containing ≥ 2 sequenced dichromatic species have enough branches to test — three do, one panel each: Trachypithecus (origin 7, 8 tips), Nomascus (origin 8, 3 tips) and Eulemur (origin 14, 2 tips); the heading of each panel gives its origin ID and the number of sequenced species (tips) it contains. The other ~12 origins are single species and cannot be tested this way.
How to read a bar. Each bar is one gene that shifted significantly in selection at that origin (Benjamini–Hochberg p_BH < 0.05). Bar length = log₂ K: bars pointing right = selection intensified, bars pointing left = selection relaxed (the ← relaxed / intensified → guides sit under each axis). Color = which module the gene belongs to (orange = pigmentation, blue = sex-hormone). Solid bar = the shift is also significant in the independent pooled analysis (corroborated); hatched bar = significant only within that one origin, so weaker evidence. A panel reading “no gene passes p(BH) < 0.05” (Eulemur) means no gene reached significance there — not that the origin is absent.
What it shows. The three origins use different, non-overlapping gene sets (quantified in §6): Trachypithecus recruits many genes from both modules, Nomascus a small pigmentation set (HRAS, POMC, HGF), and Eulemur none. There is no shared “dichromatism gene” — each origin reached the same phenotype through different parts of the coupled system.
QC. Per-origin RELAX can return an extreme boundary K on a single origin’s few branches. We cross-check every per-origin hit against the pooled RELAX fit and drop any whose extreme K is null when pooled — HPS4 is the clear case (per-origin K ≈ 30 in Trachypithecus but pooled K = 1.05, p_BH = 1.0), so it is removed as a boundary artifact. HRAS is the opposite: its per-origin K is boundary-inflated, but pooled RELAX gives a clean K = 4.3 (p_BH = 4×10⁻⁶), so it is retained and plotted at the pooled value. Bars are drawn at the pooled K wherever a per-origin K is boundary-inflated, so the figure shows de-inflated, corroborated effect sizes.
5 — Module balance per origin, corrected for panel composition
When dichromatism evolved at a given origin, was the selection concentrated in pigmentation genes, in sex-hormone genes, or split between the two? This section answers that with a single number per origin — the module balance — and shows it in the table below.
The metric. For each origin, module balance = (nP − nH) / (nP + nH), where nP and nH are the pigmentation and hormone genes under episodic selection along the origin’s branches (−1 = every selected gene is a hormone gene, +1 = every selected gene is a pigmentation gene, 0 = an even split).
The correction. The raw panel is not module-balanced by count: the hormone module was assembled as a whole endocrine pathway while the pigmentation module began as the canonical melanogenesis core, so a neutral lineage does not sit at 0. The clean-30 pigmentation expansion (NB14) brings the panel to 57 pigmentation + 53 hormone genes — a count balance of +0.036 — but for per-origin balance the correct denominator is a per-gene selection rate, not a raw count, so that an origin is not scored pigmentation-leaning simply because more pigmentation genes were tested. This section computes balance as a rate ratio; NB14’s +0.036 is the panel’s count balance and is not the same quantity (the two must not be conflated).
Table 3. Module balance for each testable origin. sigP/sigH = number of pigmentation / hormone genes under selection; rateP/rateH = those counts divided by how many genes of each module were actually tested (the panel-composition correction); rate_balance = the headline number, running from −1 (all hormone) through 0 (even) to +1 (all pigmentation). count_balance is the same idea without the correction, shown alongside so the effect of correcting is visible. NaN means no gene of either module reached significance at that origin.
Show code
# Module balance as a per-gene selection RATE per origin, using the module column the collector# adds to perorigin_v1. rate_M = (genes of module M with p_BH<0.05) / (genes of module M tested);# balance = (rate_P - rate_H) / (rate_P + rate_H), in [-1 (hormone), +1 (pigmentation)].if PER_ORIGIN isnotNoneand"module"in PER_ORIGIN.columns and"p_BH"in PER_ORIGIN.columns: bal_rows = []for o insorted(PER_ORIGIN.origin_id.unique()): sub = PER_ORIGIN[PER_ORIGIN.origin_id == o] tP = sub[sub.module =="pigmentation"].gene.nunique() tH = sub[sub.module =="hormone"].gene.nunique() sP = sub[(sub.module =="pigmentation") & (sub.p_BH <0.05)].gene.nunique() sH = sub[(sub.module =="hormone") & (sub.p_BH <0.05)].gene.nunique() rP = sP / tP if tP else0.0 rH = sH / tH if tH else0.0 bal = (rP - rH) / (rP + rH) if (rP + rH) >0elsefloat("nan") cbal = (sP - sH) / (sP + sH) if (sP + sH) >0elsefloat("nan") bal_rows.append(dict(origin=o, sigP=sP, sigH=sH, rateP=round(rP, 3), rateH=round(rH, 3), rate_balance=round(bal, 3), count_balance=round(cbal, 3))) BAL = pd.DataFrame(bal_rows)print(BAL.to_string(index=False))print("\nRate balance: +1 = purely pigmentation, -1 = purely hormone, NaN = no gene significant.")print("count_balance shown alongside to expose where the panel-composition correction matters.")print("NB14's +0.036 is the whole-panel COUNT balance (a different quantity - do not conflate).")else:print("Awaiting per-origin tables with module + p_BH columns.")
origin sigP sigH rateP rateH rate_balance count_balance
origin_14 0 0 0.000 0.00 NaN NaN
origin_7 7 6 0.135 0.12 0.057 0.077
origin_8 3 0 0.064 0.00 1.000 1.000
Rate balance: +1 = purely pigmentation, -1 = purely hormone, NaN = no gene significant.
count_balance shown alongside to expose where the panel-composition correction matters.
NB14's +0.036 is the whole-panel COUNT balance (a different quantity - do not conflate).
Figure 5 — Module balance across every origin (branch-level view)
Table 3 scores module balance only at the three RELAX-powered origins. The two figures below take the complementary branch-level view: they use the aBSREL per-tip episodic-selection calls (which exist for every dichromatic lineage, including the single-species origins RELAX cannot test), and for each lineage count how many pigmentation vs sex-hormone genes are under selection. Full captions follow each figure.
Show code
from IPython.display import Image, display_mb = os.path.join(SYN, "figures", "nb15_fig_module_balance.png")if os.path.exists(_mb): display(Image(filename=_mb))else:print("nb15_fig_module_balance.png not built - run analysis/module_selection/figures.py")
Figure 5. Each dichromatism origin tilts toward a different module.
What the data are. aBSREL episodic-selection calls (corrected p < 0.05) on the tip branches of each dichromatic origin, from results/perorigin_v1/branch_rates.csv (107-gene pigmentation + sex-hormone panel; module assigned from the panel’s own module column, not set). Built by analysis/module_selection/figures.py. Origins are grouped by family (bold headers) and, within a family, ordered by their balance.
How to read it. Each bar is one origin. Bar direction and colour give the module balance = (nP − nH)/(nP + nH): bars to the right (orange) = pigmentation-tilted, to the left (blue) = hormone-tilted, length = how lopsided. The NP/NH label on each bar is the raw count of pigmentation and hormone genes under selection at that origin. A black outline marks the three origins with ≥ 2 sequenced dichromatic species (the RELAX-powered origins of Table 3 and Fig 4); all others are single-species origins seen only at the branch level. The single-species origins sit near ±1 because they rest on very few selected genes — read those as suggestive, not as resolved module assignments.
What it shows. There is no shared module: some origins are pigmentation-tilted, others hormone-tilted, spread across the primate tree. Even the two best-powered origins differ — Trachypithecus (origin 7) leans hormone, Nomascus (origin 8) leans pigmentation — the per-origin divergence quantified in §6.
Show code
_ct = os.path.join(SYN, "figures", "nb15_fig_circular_tree_balance.png")if os.path.exists(_ct): display(Image(filename=_ct))else:print("nb15_fig_circular_tree_balance.png not built - run analysis/module_selection/tree_figures.R")
Figure 5b′. The same module balance, painted on the primate tree.
What the data are. The identical per-tip module-balance metric as Figure 5, drawn on the 117-genome primate phylogeny (circular fan layout). Built by analysis/module_selection/tree_figures.R from the same branch_rates.csv.
How to read it. Each tip carries a dot. Dot colour is the module balance on that lineage — blue = hormone-tilted, orange = pigmentation-tilted, pale = even (colour bar at centre). Dot size scales with how many genes are under selection on that lineage (few → many). Tip labels in red bold are the dichromatic species; grey labels are monochromatic. This places the balance in phylogenetic context: dichromatic lineages do not share a colour, and the module tilt varies even among close relatives (e.g. within the Trachypithecus langurs), reinforcing that each origin remodels a different part of the coupled network.
The module-balance metric (§5) scores each gene into exactly one module, and the panel assigns POMC to pigmentation (receptor_signaling; OMIM hypopigmentation phenotype). That is defensible — POMC is the precursor of α-MSH, the MC1R ligand that drives eumelanin — but it is only half of POMC’s biology. The same pro-hormone is cleaved into ACTH and β-endorphin, the HPA-axis and opioid peptides. POMC is therefore a genuine interface gene: a single locus whose products act in both the pigmentation and the endocrine modules. A one-module label is an accounting choice, not a statement that POMC is “not hormonal.”
So POMC is worth a dedicated cross-primate view, independent of which module bucket it lands in. Below we ask where in the order POMC itself is under selection — not just within the three powered origins, but across every branch the full-panel aBSREL scan covers.
What the data show. POMC shows a significant per-origin RELAX intensification in Nomascus (origin 8; K = 3.4, p_BH < 0.001 — selection intensified, not relaxed) and episodic diversifying selection (aBSREL, corrected p < 0.05) on five branches. Figure 5b shows these explicitly, colored by dichromatism state, because the pattern is mixed and that mix is the point: two of the selected tips — Nomascus concolor and N. gabriellae — are dichromatic and fall exactly on the origin-8 lineage, but a third, Macaca mulatta, is monochromatic, and the two remaining selected branches are internal Macaca clades — resolved from the POMC gene tree in the aBSREL output (one an 8-species macaque clade that is 12% dichromatic, one a 2-species clade with no dichromatic members), i.e. an essentially monochromatic macaque radiation. (HYPHY labels internal branches “NodeNN” in each gene’s own tree; those numbers are not comparable across genes, so the figure names the clade each subtends rather than the raw label.) So POMC is under selection on the origin-8 dichromatic lineage, but selection on POMC is not confined to dichromatic lineages — the same conclusion the whole-panel aBSREL scan reached (episodic selection hits dichromatic and monochromatic tips at indistinguishable rates). The reading is: POMC is an interface gene evolving in several primate lineages, one of which (the Nomascus gibbons) is a dichromatism origin. That co-localisation is a concrete follow-up target, not evidence that POMC causes the phenotype — aBSREL marks lineage-specific selection, not causation, and the Macaca signal shows the two do not track each other one-to-one.
Figure 5b — POMC selection across the primate order
The figure below asks where in the primate order POMC — a gene at the pigmentation–hormone interface — is itself under selection, and whether those branches are the dichromatic ones. Full caption follows the figure.
Show code
from IPython.display import Image, display_pomc = os.path.join(SYN, "figures", "nb15_pomc_tree.png")if os.path.exists(_pomc): display(Image(filename=_pomc))else:print("nb15_pomc_tree.png not built yet - run nb15_pomc.R")
Figure 5b. Selection on POMC across the primate order.
What the data are. POMC coding sequence from the 117 sequenced primate genomes, tested for selection with aBSREL (HyPhy), which scans every branch of the tree and flags the ones showing episodic positive selection — a burst of adaptive change confined to that branch. Results are from results/full_panel_117/absrel/POMC.ABSREL.json. POMC is singled out because it sits at the pigmentation–hormone interface (α-MSH drives pigment; ACTH/β-endorphin are hormones), so if it were the gene coupling the two modules to produce dichromatism, its selected branches should be the dichromatic ones.
How to read the two panels.(A) The 224-species tree (as in Fig 3). Each tip is colored by its coded dichromatism state (red = dichromatic, gray = monochromatic); every branch that aBSREL flags as under selection carries a star — red star if that branch is dichromatic, dark star if monochromatic. (B) The five selected branches ranked by significance; the x-axis is −log₁₀(corrected p-value), so longer bars = stronger evidence, and each bar is labeled by the branch’s dichromatism state.
What it shows. Of the five POMC-selected branches, only two are dichromatic (the Nomascus gibbons N. concolor and N. gabriellae); the other three are monochromatic macaque lineages — the M. mulatta tip and two internal Macaca clades (labeled by the clade each contains, since aBSREL’s internal “NodeNN” names are per-gene and not comparable across genes). So POMC selection overlaps one dichromatism origin (the gibbons) but is not confined to dichromatic lineages — the same lesson as the whole-panel branch scan.
Table 4. POMC’s per-origin RELAX result, read from results/perorigin_v1/per_origin_K.csv. One row per testable origin: module = the panel’s module assignment for POMC (pigmentation); K = the RELAX selection-intensity statistic (K > 1 intensified, K < 1 relaxed); p_BH = Benjamini–Hochberg-corrected p-value. POMC is significantly intensified at the Nomascus origin (origin 8, K ≈ 3.4, p_BH < 0.001) and not at the others — the per-origin counterpart to the branch-level scan in Fig 5b.
6 — Convergence and divergence across independent origins
The question. If dichromatism arose ~15–19 times, do independent origins reach it through the same genes (molecular convergence) or different ones (divergence)? This is the natural test of the two-module hypothesis: convergence would say there is one genetic route; divergence would say the coupled system can be perturbed at many points to the same phenotypic end.
Method and framing. For each powered origin the gene set is those panel genes with a significant per-origin selection-intensity shift (RELAX K, p_BH < 0.05); overlap between origins is the shared-gene count. On the current panel (size printed by the cell below; the clean-30 expansion is still landing), only two of the three powered origins carry a detectable signal — Trachypithecus (a multi-gene, both-module set) and Nomascus (a small all-pigmentation set: HGF, HRAS and POMC, the last classified as pigmentation in the panel) — and they share zero genes; Eulemur has none passing threshold. So where the comparison is well-powered it shows divergence, not convergence: independent origins are built through different genes.
This is stated as divergence-where-powered, not as a clade-wide “no convergence” claim, because the test is power-limited: with only two signalled origins, and the ~13–17 single-tip origins recoverable only by branch scans (not per-origin RELAX), a broader convergence null cannot be tested here. The overlap is recomputed live below and refreshes to 110 genes when the clean-30 expansion lands (which can only add genes to an origin’s set, so the zero-overlap divergence result is robust unless a new gene happens to hit both origins).
Table 5. Gene-set overlap between the independent origins, computed from the significant genes (p_BH < 0.05) in results/perorigin_v1/per_origin_K.csv. For each powered origin the printout lists its selected gene set, then reports the intersection between every pair of signalled origins. “ZERO overlap” means the two origins share no gene under selection — direct evidence of divergence (different molecular routes to the same phenotype) rather than convergence. Only two of the three tested origins carry a multi-gene signal, so this comparison is power-limited (see the prose above); it is reported as divergence-where-powered, not a clade-wide claim.
Show code
# Per-origin gene-set overlap from the CURRENT perorigin_v1 tables (pre-expansion run; the# printed panel size is the number of genes with RELAX results, refreshing to 110 when the# clean-30 expansion lands). Gene set per origin = genes with p_BH<0.05.import itertoolsif PER_ORIGIN isnotNoneand"p_BH"in PER_ORIGIN.columns: sig = PER_ORIGIN[PER_ORIGIN["p_BH"] <0.05] sets = {o: set(sig[sig.origin_id == o].gene) for o insorted(sig.origin_id.unique())} powered =sorted(PER_ORIGIN.origin_id.unique())print(f"powered origins tested: {powered} (>=2 dichromatic tips each)")for o in powered: g =sorted(sets.get(o, set()))print(f" {o}: {len(g)} genes under selection (p_BH<0.05) - {g if g else'(none pass threshold)'}")print() signalled = [o for o in powered if sets.get(o)]iflen(signalled) >=2:for a, b in itertools.combinations(signalled, 2): inter = sets[a] & sets[b]print(f" {a} n {b}: {len(inter)} shared genes - {sorted(inter) if inter else'ZERO overlap'}")else:print(" Only one origin carries a detectable multi-gene signal; a molecular-convergence")print(" test needs >=2 signalled origins, so the comparison is power-limited (see prose).")print(f"\n(Current: {PER_ORIGIN['gene'].nunique()}-gene panel. Refreshes when clean-30 lands -> 110 genes.)")else:print("Awaiting per-origin tables with p_BH for the overlap recompute.")
Figure 6 — The dichromat-vs-sister contrast, gene by gene
Table 5 counts shared genes between origins. The two figures below make the underlying contrast explicit at the finest resolution the data allow: for each dichromatic taxon, they compare its selected genes against those of its closest monochromatic relative on the tree. If a gene is selected in the dichromat but not its sister, that lineage-specific difference is a candidate for where the phenotype’s genetic basis sits (a difference marks lineage-specific selection, not proven causation). Full captions follow each figure.
Show code
_sc = os.path.join(SYN, "figures", "nb15_fig_sister_pair_contrast.png")if os.path.exists(_sc): display(Image(filename=_sc))else:print("nb15_fig_sister_pair_contrast.png not built - run analysis/module_selection/sister_figures.py")
Figure 6. Selection in each dichromatic taxon vs its closest monochromatic sister, by gene.
What the data are. aBSREL episodic-selection calls (corrected p < 0.05) from results/perorigin_v1/branch_rates.csv (107-gene panel), for 21 dichromatic species and their nearest monochromatic sister taxon (sister pairs in analysis/module_selection/sister_pairs.csv, taken from the phenotype tree). Built by analysis/module_selection/sister_figures.py.
How to read it. The figure is split into two stacked panels so every gene label is legible: the pigmentation module (top, orange) and the sex-hormone module (bottom, blue). Both panels share the same rows — dichromatic taxa grouped by clade (bold clade headers; the small grey line under each taxon names its monochromatic sister). Columns are the genes of that module. Each cell is split on the diagonal: the lower-left triangle = the dichromatic species, upper-right triangle = its sister. A solid triangle = that gene is under episodic selection in that taxon; a faded/pale triangle = selected in the sister only; an empty grey cell = selected in neither. So a solid-orange lower-left triangle with an empty upper-right (top panel) is a pigmentation gene selected in the dichromat but not its sister — a lineage-specific difference.
Margins (added for counting). In each panel the right margin gives that module’s per-pair totals — the number of genes under selection in the dichromat (D) and in its sister (M) for that module; the pair is printed in red where the dichromat has more. The bottom margin gives each gene’s column totals — the number of sister-pairs in which that gene is selected in the dichromat (red, D) vs the sister (grey, M). The italic line under the title sums the two modules across all 21 pairs.
What it shows. The filled cells do not line up into shared columns: different dichromatic taxa light up different genes, and few genes are selected in a dichromat and absent in its sister across more than one clade. The Trachypithecus langurs (a single origin, many sequenced species) form the densest block and share genes among themselves, but that block does not repeat in the gibbons, guenons, howlers or lemurs — the by-eye version of the zero cross-origin overlap in Table 5.
Counting the contrast: is there more selection in dichromats or in their sisters?
The figure’s margins invite a direct count. Summed over all 21 sister-pairs, and split by module, the totals answer whether episodic selection is concentrated on the dichromatic lineages.
Show code
# Reproduce the exact per-cell selection calls the figure uses, then total them by module._SP = load_input("sister_pairs")_GM = load_input("gene_modules")_gmod =dict(zip(_GM.gene, _GM.module))def _selset(sp): # genes under episodic selection on that tip branch (aBSREL corrected p<0.05)returnset(BRANCH[(BRANCH.branch == sp) & (BRANCH.selected_flag ==True) & (BRANCH.is_tip ==True)].gene)_rows = []for _, _r in _SP.iterrows(): _sd, _sm = _selset(_r.dich), _selset(_r.sister_mono)def _mc(gs, mod): returnsum(1for g in gs if _gmod.get(g) == mod) _rows.append(dict(dichromat=_r.dich.replace("_", " "), sister=_r.sister_mono.replace("_", " "), dich_pig=_mc(_sd, "pigmentation"), dich_hor=_mc(_sd, "hormone"), mono_pig=_mc(_sm, "pigmentation"), mono_hor=_mc(_sm, "hormone")))_cnt = pd.DataFrame(_rows)_summary = pd.DataFrame({"module": ["pigmentation", "hormone", "ALL"],"dichromat_selection_events": [_cnt.dich_pig.sum(), _cnt.dich_hor.sum(), _cnt.dich_pig.sum() + _cnt.dich_hor.sum()],"monochromatic_sister_events": [_cnt.mono_pig.sum(), _cnt.mono_hor.sum(), _cnt.mono_pig.sum() + _cnt.mono_hor.sum()],})_summary["difference_(dich-mono)"] = (_summary.dichromat_selection_events- _summary.monochromatic_sister_events)print("Selection events across", len(_SP), "dichromatic-vs-sister pairs:\n")print(_summary.to_string(index=False))
[sister_pairs] comparative-genomics/analysis/module_selection/sister_pairs.csv
what : Each dichromatic taxon paired with its closest monochromatic relative on the tree
source : Nearest-monochromatic-sister search on the phenotype tree
produced : analysis/module_selection (sister-pair construction)
loaded : 21 rows x 5 cols | sha256[:16]=a7ebce55618c4c9f
[gene_modules] comparative-genomics/analysis/module_selection/data/gene_modules.csv
what : Module label (pigmentation / hormone) for every gene in the 107-gene panel
source : Panel definition (NB13/NB14)
produced : analysis/module_selection (panel module map)
loaded : 107 rows x 2 cols | sha256[:16]=6495eb55dd747735
Selection events across 21 dichromatic-vs-sister pairs:
module dichromat_selection_events monochromatic_sister_events difference_(dich-mono)
pigmentation 61 72 -11
hormone 51 49 2
ALL 112 121 -9
Table 5b. Episodic-selection events in dichromatic taxa vs their monochromatic sisters, by module.
What the data are. Counts of gene-level episodic-selection calls (aBSREL corrected p < 0.05, results/perorigin_v1/branch_rates.csv) summed over the 21 dichromatic-vs-sister pairs in analysis/module_selection/sister_pairs.csv, split by gene module (gene_modules.csv). A “selection event” is one gene under selection on one taxon’s terminal branch; the same gene counts once per taxon it is selected in.
What it shows. Selection is not concentrated on the dichromatic lineages: dichromats carry 112 selection events (61 pigmentation + 51 hormone) against 121 in their monochromatic sisters (72 pigmentation + 49 hormone). The dichromats have marginally fewer pigmentation events and a near-identical hormone count. This is the pairwise, module-resolved counterpart to the phylogenetically-controlled test (whole-tree contrast: null after correcting for tree structure) — the difference that makes a lineage dichromatic is in which genes are selected (Figure 6), not in how many.
Show code
_sn = os.path.join(SYN, "figures", "nb15_fig_sister_network_diff.png")if os.path.exists(_sn): display(Image(filename=_sn))else:print("nb15_fig_sister_network_diff.png not built - run analysis/module_selection/sister_figures.py")
Figure 6b′. Where the dichromat-vs-sister difference sits in the coupled network.
What the data are. The same dichromat-vs-sister selection contrast as Figure 6, for four representative pairs, drawn on the coupled pigmentation–hormone gene network instead of a matrix. Node positions come from a fixed layout (network_layout.json) so genes sit in the same place in every panel; edges are human protein–protein interactions from STRING v12 (score ≥ 0.4). Built by analysis/module_selection/sister_figures.py.
How to read it. Each of the four panels is one dichromatic taxon (named top-left) vs its sister. Genes are nodes: the pigmentation module fills the left (cream) box, the hormone module the right (pale-blue) box. A node is red-ringed and filled = selected in the dichromatic species only (the lineage-specific variation), grey-filled = selected in both taxa, open/hollow = selected in the sister only, and faint = selected in neither. Because the layout is shared, you can scan the same gene across panels.
What it shows. The red-ringed (dichromat-specific) genes fall in different parts of the network in each pair — Trachypithecus francoisi lights up a broad pigmentation+hormone spread, Nomascus concolor a pigmentation-side cluster around POMC/SOX10, Cercopithecus hamlyni a small pigmentation set, Alouatta caraya a hormone-side cluster — again the divergence story: the same phenotype, reached by perturbing different nodes of one coupled system.
7 — Synthesis
Sexual dichromatism in primates is a labile, polygenic trait built from two coupled modules. It has arisen on the order of 15–19 times and is lost far faster than it is gained (§3). At the two origins where its genetic basis can actually be tested, the signalled gene sets do not overlap (§4, §6, Fig 6) and differ in module balance (§5, Fig 5) — divergence where the comparison is powered. Whether origins differ systematically in architecture across the clade is an open question the current sampling cannot resolve (only two of ~15–19 origins carry a per-origin signal), so it is framed as underpowered rather than claimed.
What the sequence data add, beyond the tree. The lability and origin count (§3) come from phenotype and phylogeny alone. The genetics answers what those cannot: which molecular change accompanied each origin. Three results here require sequence data and would be invisible to a tree-only analysis — (1) MC1R, the gene birds use again and again for dichromatism, is not the primate route (no pooled signal, K = 0.45; it surfaces in one lineage only); (2) the two testable origins recruited non-overlapping gene sets from both the pigmentation and the hormone module (§4, §6); and (3) POMC, a gene that sits at the pigment–hormone interface, is itself under selection at one origin (§5b). Together these say a trait built one way in birds is built many different ways in primates, drawing on a coupled pigmentation–hormone system rather than a single master gene.
This is exactly what a network framing predicts and a single-gene framing would miss: if dichromatism is the output of a coupled system, there are many points at which it can be pushed to produce (or lose) the phenotype, and different lineages have used different ones. The panel that makes this testable — which genes, and why — is built in NB14, with the hormone layer developed in NB13. The selection evidence is strongest for the three powered origins; the single-tip origins and the convergence test are power-limited and framed as such.
