Module 6: Phylogeny & Evolution of Flowers

Darwin called the sudden Cretaceous appearance of angiosperms “an abominable mystery.” This module reconstructs flower evolution: the deep-branching ANA grade, the explosive radiation of eudicots, molecular-clock dating, convergent pollination syndromes, and the extraordinary diversification of the orchid family — 28,000 species from a single Cretaceous ancestor.

1. Origin of Flowers — Darwin’s “Abominable Mystery”

In an 1879 letter to Joseph Hooker, Darwin complained that “the rapid development, as far as we can judge, of all the higher plants within recent geological times is an abominable mystery.” The fossil record shows the nearly simultaneous appearance of diverse angiosperm lineages in the Early Cretaceous (~140 Ma), with no clear extinct stem group.

Today we recognise unambiguous angiosperm pollen and leaves from the Hauterivian (~130 Ma) and tentative records back to the Valanginian or even Late Jurassic. Controversial Triassic pollen ( Crinopolles, Hochuli & Feist-Burkhardt 2013 ) could push the origin to ~200 Ma, but most morphological and molecular evidence supports a Jurassic stem lineage and a Cretaceous crown radiation.

Why flowers won

Flowering plants replaced gymnosperms across most of the Earth’s surface in under 80 million years. Candidate drivers include:

Double fertilization & endosperm — efficient resource packaging for offspring.
Animal pollination — directional, long-distance gene flow and rapid isolation.
Fruit dispersal — co-evolution with vertebrates expanded seed shadows.
Vein density and leaf hydraulics — higher transpiration, faster carbon gain (Boyce et al. 2009).
Genome duplications — whole-genome duplication events immediately precede major radiations.

2. Angiosperm Phylogeny & the ANA Grade

Molecular phylogenies (Soltis et al. 1999; APG IV 2016) resolved the root of the angiosperm tree with three successive sister lineages — collectively the ANA grade:

Amborella trichopoda — a single shrub species endemic to New Caledonia, sister to all other living angiosperms.
Nymphaeales — water lilies, including Nymphaea and Cabomba.
Austrobaileyales — including Illicium (star anise) and Trimenia.

The rest of the tree splits into Monocots (~60,000 species) and Eudicots (~175,000 species). The magnoliids are a smaller, basal clade sister to either monocots or eudicots depending on the analysis.

The eudicot explosion

Triaperturate (three-pore) pollen — the eudicot synapomorphy — first appears in the Barremian (~125 Ma). Within 20 Myr, the core eudicot split into rosids and asterids, which together comprise ~75% of modern angiosperm species.

Derivation: molecular clock divergence times

Under the strict molecular clock, the expected number of substitutions per site between two lineages that diverged at time \(t\) and each evolve at rate\(\mu\) is

\( K = 2\,\mu\,t \quad\Longrightarrow\quad t = \dfrac{K}{2\mu} \)

For plastid \(\text{rbcL}\), \(\mu \approx 3.2\times10^{-10}\) subs/site/yr.

Real lineages violate the strict clock, so relaxed-clock Bayesian methods (Thorne & Kishino 2002; BEAST) place a prior on per-lineage rates, typically log-normal:

\( \log \mu_b \sim \mathcal{N}(\log \bar{\mu},\,\sigma^2) \)

Fossil calibrations convert the tree from substitutions to years.

Applying this to a matrix of chloroplast genes across the ANA grade and the major clades yields crown ages for angiosperms clustered around 160–180 Ma, with 95% credible intervals extending to 210 Ma (Magallón et al. 2015).

Angiosperm Phylogeny with Key Innovations

3. Convergent Evolution of Pollination Syndromes

The same suite of flower traits — long red corolla tubes, scentlessness, diurnal opening, copious dilute nectar — has evolved independently in scores of unrelated lineages pollinated by hummingbirds. This is a textbook case of convergent evolution.

Testing convergence statistically

Given a phylogenetic tree and a binary trait, the minimum number of trait origins is found by Fitch’s parsimony algorithm. Under a neutral model where the trait evolves with rate \(q\) (symmetric Markov), the expected number of changes along a tree of total branch length \(L\) is \(qL\). Observed numbers substantially larger than permuted nulls indicate selection for the trait state.

\( P(\text{origins} \ge k \mid \text{null}) \;=\; \dfrac{1}{B}\sum_{b=1}^{B} \mathbb{1}[\text{origins}(\pi_b) \ge k] \)

Permutation test: shuffle trait labels \(B\) times, count origins each time.

Iconic convergences

Hummingbird syndrome: Costus, Heliconia, Lobelia, Penstemon, Columnea, Ipomopsis — 30+ independent origins in the Americas alone.
Bat syndrome: Agave, Musa, Ceiba, many cacti — pale, nocturnal, fermented-smelling flowers with abundant nectar.
Hawkmoth syndrome: Oenothera, Nicotiana, Datura, Angraecum — long white tubes, evening-scented.
Carrion syndrome: Rafflesia, Stapelia, Amorphophallus — red/purple, putrid, fly-pollinated.

4. The Orchid Radiation

Orchidaceae is the largest angiosperm family by species count: ~28,000 species in ~750 genera, or roughly one in ten of all flowering plants. The crown age of the family is only ~80 Myr, implying an extraordinary net diversification rate.

Derivation: net diversification rate

Under a pure-birth (Yule) model, the number of lineages at time \(t\) is\(N(t) = N_0\,e^{r t}\). For Orchidaceae, \(N_0 = 1\) at crown age \(t \approx 80\) Myr and \(N = 28{,}000\), giving

\( r = \dfrac{\ln(N/N_0)}{t} = \dfrac{\ln 28000}{80} \approx 0.13\;\text{Myr}^{-1} \)

About 2–3× the background angiosperm rate (~0.05 Myr^-1).

What fuelled the orchid radiation?

Epiphytism (evolved ~40 Ma in Epidendroideae) opened a three-dimensional canopy niche.
Pollinia package thousands of pollen grains for single-visit delivery, enabling extreme pollinator specialization.
Deceptive pollination (food deception, sexual deception, brood-site deception) reduces the need to reward visitors and favours rapid floral divergence.
Mycoheterotrophic seedlings: dust-like seeds with no food reserves, germinating only after symbiosis with a compatible fungus.
Resupination: twisting of the pedicel through 180° produces the lip-down architecture that defines orchid flowers.

Givnish et al. (2015) showed that each of these innovations is associated with a step-change in diversification rate, and epiphytic orchid clades have rates up to\(r \approx 0.18\) Myr^-1, among the highest in angiosperms.

Simulation: Molecular Clock Divergence Times

Relaxed log-normal clock applied to ANA-grade splits, lineage-through-time plot across 140 Myr, rate heterogeneity between plastid, nuclear, and mitochondrial markers, and the Bayesian posterior on the crown angiosperm age.

Python

script.py117 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Molecular clock divergence time estimation for angiosperm lineages.
# Simple strict clock: t = K / (2 * mu)   where K = observed sequence divergence
# Relaxed log-normal clock: per-lineage rate ~ LogNormal(mu, sigma).
# Calibrations: oldest angiosperm-like pollen ~ 140 Ma (Hauterivian).

np.random.seed(1)

# Observed pairwise divergences (substitutions per site, chloroplast rbcL)
# Synthetic but calibrated to typical published values (e.g. Magallon & Castillo 2009).
pairs = {
    'Amborella vs core angios': 0.113,
    'Nymphaeales vs core': 0.102,
    'Austrobaileyales vs core': 0.095,
    'Monocots vs eudicots': 0.082,
    'Asterid vs rosid': 0.064,
    'Orchidaceae crown': 0.032,
    'Asteraceae crown': 0.028,
}

mu_rbcL = 3.2e-10  # substitutions per site per year (rbcL)
times_strict = {k: v / (2 * mu_rbcL) / 1e6 for k, v in pairs.items()}

# Relaxed clock: draw per-lineage rates log-normally, take mean divergence time
sigma = 0.35
n_draws = 2000
times_relaxed = {}
times_ci = {}
for name, K in pairs.items():
    rates = mu_rbcL * np.exp(np.random.normal(0, sigma, n_draws))
    t_draws = K / (2 * rates) / 1e6
    times_relaxed[name] = np.median(t_draws)
    times_ci[name] = (np.percentile(t_draws, 2.5), np.percentile(t_draws, 97.5))

# Angiosperm diversification pattern (Magallon et al. 2015 style)
# Lineage-through-time: N(t) = N0 * exp(r*t) with episodic rate shifts
t_grid = np.linspace(0, 140, 500)  # Ma BP -> present if we reverse
# Simple pure-birth with rate 0.093 /Myr after Cretaceous radiation
lambda_bg = 0.008
lambda_radiation = 0.093
N_lineages = np.ones_like(t_grid)
for i in range(1, len(t_grid)):
    dt = t_grid[i] - t_grid[i - 1]
    lam = lambda_radiation if t_grid[i] > 20 and t_grid[i] < 110 else lambda_bg
    N_lineages[i] = N_lineages[i - 1] * np.exp(lam * dt)

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0a0a1a')
for ax in (ax1, ax2, ax3, ax4):
    ax.set_facecolor('#0a0a1a')
    ax.tick_params(colors='#94a3b8')
    for sp in ax.spines.values():
        sp.set_color('#334155')
    ax.grid(True, color='#334155', alpha=0.35)

names = list(times_strict.keys())
ys = np.arange(len(names))[::-1]
med = [times_relaxed[n] for n in names]
lo = [times_ci[n][0] for n in names]
hi = [times_ci[n][1] for n in names]

ax1.errorbar(med, ys, xerr=[np.array(med) - np.array(lo), np.array(hi) - np.array(med)],
             fmt='o', color='#f472b6', ecolor='#fb7185',
             markersize=10, linewidth=2, capsize=5)
ax1.set_yticks(ys)
ax1.set_yticklabels(names, color='#cbd5e1')
ax1.set_xlabel('Divergence time (Ma)', color='#cbd5e1')
ax1.set_title('A. Relaxed molecular clock estimates (95% CI)',
              color='#f9a8d4', fontweight='bold')
ax1.invert_xaxis()

ax2.plot(t_grid, N_lineages, color='#f472b6', linewidth=2.5)
ax2.fill_between(t_grid, N_lineages, alpha=0.2, color='#f472b6')
ax2.axvspan(110, 20, alpha=0.1, color='#fde68a', label='Cretaceous-Paleogene radiation')
ax2.set_yscale('log')
ax2.set_xlabel('Time (Ma, present on left)', color='#cbd5e1')
ax2.set_ylabel('log(lineages)', color='#cbd5e1')
ax2.set_title('B. Lineage-through-time plot (angiosperms)',
              color='#f9a8d4', fontweight='bold')
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#e2e8f0')
ax2.invert_xaxis()

# Rate heterogeneity across rbcL / ITS / mitochondrial
markers = ['rbcL (cpDNA)', 'ITS (nuclear)', 'matR (mtDNA)']
rates = [3.2e-10, 2.8e-9, 1.2e-10]
ax3.bar(markers, np.array(rates) * 1e9, color=['#f472b6', '#fb7185', '#fde68a'],
        edgecolor='white', linewidth=0.5)
ax3.set_ylabel('Rate (10^-9 subs / site / yr)', color='#cbd5e1')
ax3.set_title('C. Rate heterogeneity across loci',
              color='#f9a8d4', fontweight='bold')

# Bayesian posterior on crown angiosperm age
ages = np.random.normal(168, 35, 8000)
ages = np.clip(ages, 100, 280)
ax4.hist(ages, bins=50, color='#f472b6', alpha=0.75, edgecolor='white', linewidth=0.3)
ax4.axvline(140, color='#fde68a', linestyle='--', linewidth=2,
            label='Oldest unambiguous fossil (140 Ma)')
ax4.axvline(np.median(ages), color='#fb7185', linestyle='-', linewidth=2,
            label=f'Median posterior {np.median(ages):.0f} Ma')
ax4.set_xlabel('Crown angiosperm age (Ma)', color='#cbd5e1')
ax4.set_ylabel('Posterior density', color='#cbd5e1')
ax4.set_title('D. Age of the angiosperm crown (Bayesian posterior)',
              color='#f9a8d4', fontweight='bold')
ax4.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#e2e8f0')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
print('Molecular clock divergence times:')
for name, t in times_relaxed.items():
    lo_, hi_ = times_ci[name]
    print(f'  {name}: {t:.0f} Ma (95% CI {lo_:.0f} - {hi_:.0f})')
print(f'Crown angiosperm median: {np.median(ages):.0f} Ma')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: Convergent Evolution of Hummingbird Flowers

Simulated phylogeny of 400 lineages; permutation test for excess clade-origins of the hummingbird pollination syndrome; convergent elongation of corolla tubes; Dollo gain:loss asymmetry.

Python

script.py107 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Convergent evolution of hummingbird-pollinated flowers across angiosperms.
# Simulate independent origins on a phylogenetic tree and test whether
# the observed number of transitions exceeds neutral expectation.

np.random.seed(6)

# 1. Simulated phylogeny of 400 plant lineages, each scored for:
#    - flower colour (0 = other, 1 = red/orange)
#    - corolla tube length (mm)
#    - pollinator syndrome (0 = non-hummer, 1 = hummer)

n = 400
colour = np.random.binomial(1, 0.08, n)
tube_length = np.where(colour == 1,
                       np.random.gamma(3.5, 7, n),
                       np.random.gamma(1.5, 4, n))
# Hummingbird syndrome probability given (colour, tube length)
p_hb = 1.0 / (1.0 + np.exp(-(0.6 * colour + 0.05 * tube_length - 3.5)))
syndrome = np.random.binomial(1, p_hb)

# 2. Observed origins: cluster species by a naive phylogenetic signal
# Assume a tree with 20 clades of 20 species each; count clades that
# contain >= 1 hummingbird-syndrome species.
clade = np.repeat(np.arange(20), 20)
observed_origins = 0
for c in range(20):
    if np.any(syndrome[clade == c]):
        observed_origins += 1

# 3. Null distribution: shuffle syndrome labels, re-count origins
nulls = []
for _ in range(2000):
    sh = np.random.permutation(syndrome)
    origins = sum(1 for c in range(20) if np.any(sh[clade == c]))
    nulls.append(origins)
nulls = np.array(nulls)
p_value = np.mean(nulls >= observed_origins)

# 4. Tube-length convergence: compare means by syndrome across clades
means_hb = []
means_non = []
for c in range(20):
    mask = clade == c
    if np.sum(syndrome[mask]) > 0:
        means_hb.append(np.mean(tube_length[mask & (syndrome == 1)]))
    if np.sum(syndrome[mask]) < np.sum(mask):
        means_non.append(np.mean(tube_length[mask & (syndrome == 0)]))

# 5. Rate of convergent evolution: Dollo lambda
# Fraction of trait transitions that are independent gains
gains = observed_origins
losses = max(0, int(0.25 * gains))

ax1.scatter(tube_length, colour + 0.05 * np.random.randn(n),
            c=syndrome, cmap='RdPu', s=30, alpha=0.7, edgecolor='white',
            linewidth=0.3)
ax1.set_xlabel('Corolla tube length (mm)', color='#cbd5e1')
ax1.set_ylabel('Red/orange colour (jittered)', color='#cbd5e1')
ax1.set_title('A. Trait covariation across lineages',
              color='#f9a8d4', fontweight='bold')

ax2.hist(nulls, bins=30, color='#fb7185', alpha=0.55, edgecolor='white',
         linewidth=0.3, label='null (permuted)')
ax2.axvline(observed_origins, color='#fde68a', linewidth=3,
            label=f'Observed = {observed_origins} origins (p = {p_value:.3f})')
ax2.set_xlabel('Number of clade origins', color='#cbd5e1')
ax2.set_ylabel('Count', color='#cbd5e1')
ax2.set_title('B. Null test for excess convergence',
              color='#f9a8d4', fontweight='bold')
ax2.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#e2e8f0')

ax3.boxplot([means_non, means_hb], labels=['non-hummer', 'hummer'],
            patch_artist=True,
            boxprops=dict(facecolor='#f472b6', color='#fb7185'),
            medianprops=dict(color='#fde68a'))
ax3.set_ylabel('Corolla tube length (mm)', color='#cbd5e1')
ax3.set_title('C. Convergent elongation in hummer clades',
              color='#f9a8d4', fontweight='bold')

# Dollo ratio bar
ax4.bar(['Gains', 'Losses'], [gains, losses],
        color=['#f472b6', '#fb7185'], edgecolor='white', linewidth=0.5)
ax4.set_ylabel('Events (across tree)', color='#cbd5e1')
ax4.set_title(f'D. Gain:loss ratio = {gains}/{max(losses,1)}',
              color='#f9a8d4', fontweight='bold')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
print(f'Observed clade origins of hummer syndrome: {observed_origins} / 20')
print(f'Null mean: {nulls.mean():.2f}, permutation p-value: {p_value:.4f}')
print('Red tubular flowers evolved independently many times -')
print('a classic case of convergent selection by shared pollinators.')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Simulation: Diversification Rate Analysis of Orchids

Birth-death simulations contrasting orchid-like and background diversification rates; species-per-subfamily distribution; key innovations (epiphytism, pollinia, deception); age vs. richness scatter for major angiosperm families.

Python

script.py103 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Diversification rate analysis for the orchid radiation.
# Birth-death model:  dN/dt = (lambda - mu) N = r N
# Net rate r = 0.05-0.10 /Myr for angiosperms overall, but orchids
# appear to have r ~ 0.18 /Myr on some branches.

# 1. Lineage-through-time simulation of a birth-death tree
def simulate_bd(lam, mu, T):
    t = 0.0
    lineages = [0.0]
    n = 1
    times, counts = [0.0], [1]
    while t < T and n > 0:
        rate = (lam + mu) * n
        dt = np.random.exponential(1.0 / rate)
        t += dt
        if t > T:
            break
        if np.random.random() < lam / (lam + mu):
            n += 1
        else:
            n -= 1
        times.append(t)
        counts.append(n)
    return np.array(times), np.array(counts)

np.random.seed(11)
T = 80.0
lam1, mu1 = 0.18, 0.05       # orchid-like
lam2, mu2 = 0.08, 0.05       # background angiosperm
t_orch, n_orch = simulate_bd(lam1, mu1, T)
t_bg, n_bg = simulate_bd(lam2, mu2, T)

# 2. Species richness across major orchid tribes (approximate, Chase 2005)
tribes = ['Apostasioideae', 'Vanilloideae', 'Cypripedioideae',
          'Orchidoideae', 'Epidendroideae']
species_counts = [15, 180, 170, 3600, 22000]

# 3. Key innovations and species density
traits = ['Epiphyty', 'Pollinia', 'Mycoheterotrophy',
          'Deceptive pollination', 'Resupination']
diversity_effect = [1.9, 1.2, 0.4, 1.5, 1.1]   # fold increase in richness

# 4. Age vs species richness for angiosperm families
np.random.seed(4)
ages = np.array([35, 55, 70, 85, 100, 115, 130, 130, 140])
families = ['Orchidaceae', 'Asteraceae', 'Fabaceae', 'Poaceae', 'Rubiaceae',
            'Solanaceae', 'Ranunculaceae', 'Nymphaeaceae', 'Amborellaceae']
richness = [28000, 33000, 19500, 12000, 13000, 2700, 2500, 70, 1]

ax1.step(t_orch, n_orch, color='#f472b6', linewidth=2,
         label=f'orchid-like (r={lam1 - mu1:.2f})')
ax1.step(t_bg, n_bg, color='#fb7185', linewidth=2,
         label=f'background (r={lam2 - mu2:.2f})')
ax1.set_yscale('log')
ax1.set_xlabel('Time (Myr)', color='#cbd5e1')
ax1.set_ylabel('log(lineages)', color='#cbd5e1')
ax1.set_title('A. Simulated birth-death lineage-through-time',
              color='#f9a8d4', fontweight='bold')
ax1.legend(facecolor='#1e293b', edgecolor='#334155', labelcolor='#e2e8f0')

ax2.barh(tribes, species_counts, color='#f472b6', edgecolor='white', linewidth=0.5)
ax2.set_xscale('log')
ax2.set_xlabel('Species count (log)', color='#cbd5e1')
ax2.set_title('B. Species per orchid subfamily',
              color='#f9a8d4', fontweight='bold')

ax3.bar(traits, diversity_effect, color='#fb7185', edgecolor='white', linewidth=0.5)
ax3.set_ylabel('Fold richness increase', color='#cbd5e1')
ax3.set_title('C. Key innovations and orchid diversification',
              color='#f9a8d4', fontweight='bold')
plt.setp(ax3.get_xticklabels(), rotation=20, ha='right')

ax4.scatter(ages, richness, s=120, color='#f472b6', edgecolor='white', linewidth=0.8)
for age, r, f in zip(ages, richness, families):
    ax4.annotate(f, (age, r), color='#cbd5e1', fontsize=8,
                 xytext=(4, 4), textcoords='offset points')
ax4.set_yscale('log')
ax4.set_xlabel('Crown age (Ma)', color='#cbd5e1')
ax4.set_ylabel('Species richness (log)', color='#cbd5e1')
ax4.set_title('D. Age vs richness (no simple relation)',
              color='#f9a8d4', fontweight='bold')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a1a')
print(f'Simulated orchid-like lineage count at {T:.0f} Myr: {n_orch[-1]}')
print(f'Simulated background lineage count at {T:.0f} Myr: {n_bg[-1]}')
print('Real orchid family has ~28,000 species - ~10% of all flowering plants.')
print('Age is a poor predictor of richness; key innovations dominate.')

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key References

• Darwin, C. (1879). Letter to J. D. Hooker, 22 July 1879 (“an abominable mystery”).

• APG IV (2016). “An update of the Angiosperm Phylogeny Group classification for the orders and families of flowering plants.” Bot. J. Linn. Soc., 181, 1–20.

• Soltis, P. S., Soltis, D. E. & Chase, M. W. (1999). “Angiosperm phylogeny inferred from multiple genes as a tool for comparative biology.” Nature, 402, 402–404.

• Magallón, S. et al. (2015). “A metacalibrated time-tree documents the early rise of flowering plant phylogenetic diversity.” New Phytologist, 207, 437–453.

• Boyce, C. K. et al. (2009). “Angiosperm leaf vein evolution was physiologically and environmentally transformative.” Proc. R. Soc. B, 276, 1771–1776.

• Hochuli, P. A. & Feist-Burkhardt, S. (2013). “Angiosperm-like pollen and Afropollis from the Middle Triassic.” Frontiers in Plant Science, 4, 344.

• Thorne, J. L. & Kishino, H. (2002). “Divergence time and evolutionary rate estimation with multilocus data.” Systematic Biology, 51, 689–702.

• Fenster, C. B. et al. (2004). “Pollination syndromes and floral specialization.” Annu. Rev. Ecol. Evol. Syst., 35, 375–403.

• Givnish, T. J. et al. (2015). “Orchid phylogenomics and multiple drivers of their extraordinary diversification.” Proc. R. Soc. B, 282, 20151553.

• Chase, M. W. et al. (2015). “An updated classification of Orchidaceae.” Bot. J. Linn. Soc., 177, 151–174.

• Friis, E. M., Crane, P. R. & Pedersen, K. R. (2011). Early Flowers and Angiosperm Evolution. Cambridge UP.

← Module 5: Reproduction & Genetics Module 7: Ornamental Horticulture →