What is spatial phylogenetics?

Why measure biodiversity using phylogeny?

Why measure biodiversity using phylogeny?

The unit of phylogenetic diversity is the clade

Every branch becomes a community matrix column

Add each clade to the species-level community matrix as a new “taxon”, and store its branch length as a weight. Diversity metrics then work on branches just as they do on species: PD is the sum of the weights of the branches present.

Common misconceptions

Phylogenetic endemism

Endemism weights each taxon by the inverse of its range size, so restricted taxa count more; phylogenetic endemism (PE) applies this to branches, summing each branch’s length divided by its range size.

Ranges don’t have to be binary

A clade’s value comes from its members: present if any are (binary), P = 1 − ∏(1 − p_i) (probability), or the sum (abundance).

Quantitative measures paint richer pictures

Moss neo- and paleo-endemism (Kling et al. 2026, Systematic Botany).

The phylospatial package

Other options for PD analyses in R

Many R packages compute PD across sites (Kling 2025, MEE, Table 1)

Two ingredients

phylospatial() combines a tree and community data into one object, building a range for every branch.

A toy dataset in one line

ps_simulate() generates a random dataset, useful for learning and testing.

sim <- ps_simulate(n_tips = 20, data_type = "probability")
plot(sim, "comm", taxa = 1:4, nc = 4, axes = FALSE)

From occurrence points to a grid

ps_grid() turns occurrence records into an equal-area community grid. Raw records reflect sampling effort; modeled ranges smooth it out.

pts    <- read.csv("data/moss_points.csv")  # species, lon, lat
comm   <- ps_grid(pts)                      # equal-area grid
ps_pts <- phylospatial(comm, ps$tree)

Building the moss dataset

443 California moss species: occurrence probabilities from species distribution models on a 30 km grid, plus a phylogeny (Kling et al. 2026, Systematic Botany). Both ship with the package; moss() is a one-line shortcut for the code below.

ext  <- system.file("extdata", package = "phylospatial")
comm <- rast(file.path(ext, "moss_comm.tif"))
tree <- ape::read.tree(file.path(ext, "moss_tree.nex"))
ps   <- phylospatial(comm, tree)

Community data type detected: probability

ps

`phylospatial` object
  - 884 lineages across 527 occupied sites (1116 total) 
  - community data type: probability 
  - branch length rescaling: sum1 
  - site weights: equal 
  - spatial data class: SpatRaster 
  - dissimilarity data: none 

plot(ps, "tree", type = "fan",
     show.tip.label = FALSE)

Inside a phylospatial object

str(ps, max.level = 1)

List of 10
 $ comm       : num [1:527, 1:884] 0.0154 0.0271 0.0582 0.2086 0.2531 ...
  ..- attr(*, "dimnames")=List of 2
 $ tree       :List of 4
  ..- attr(*, "class")= chr "phylo"
  ..- attr(*, "order")= chr "cladewise"
 $ spatial    :S4 class 'SpatRaster' [package "terra"]
 $ occupied   : int [1:527] 1 2 3 4 5 6 7 8 9 10 ...
 $ n_sites    : int 1116
 $ site_weight: NULL
 $ data_type  : chr "probability"
 $ clade_fun  :function (x)  
  ..- attr(*, "srcref")= 'srcref' int [1:8] 129 30 129 56 30 56 3623 3623
  .. ..- attr(*, "srcfile")=Classes 'srcfilealias', 'srcfile' <environment: 0x97aa4dab8> 
 $ rescale    : chr "sum1"
 $ dissim     : NULL
 - attr(*, "class")= chr "phylospatial"

Diversity and endemism metrics

ps_diversity() computes various measures of richness, divergence, and regularity, each in phylogenetic and species-level variations (?ps_diversity for the full list).

div <- ps_diversity(ps, metric = c("TD", "PD", "PE"))

Richness: species richness, PD, and PE

Each map scaled to its maximum. PD closely tracks richness; PE shifts weight toward places where range-restricted lineages concentrate.

Divergence and regularity

Beyond how much evolutionary history a site holds, ps_diversity() measures how it’s structured. Divergence (MPD): how distantly related co-occurring lineages are, on average. Regularity (VPD): how evenly spaced they are across the tree; low variance means evenly spaced lineages, high means clumped clusters of close relatives.

struct <- ps_diversity(ps, metric = c("MPDT", "VPDT"))

Which part of evolutionary history?

Recall that branch lengths can measure divergence, time, or diversification events. Tree scaling functions change which of those, and which part of history, the metrics weigh. Put the result in ps$tree and every metric follows.

tr <- rescale_tree(ps$tree, "tip1")   # root-to-tip height = 1
slice_tree(tr, min = 0.5)             # older half only
slice_tree(tr, max = 0.5)             # recent half only
delta_tree(tr, delta = 0.3)           # emphasize deep branches
delta_tree(tr, delta = 3)             # emphasize recent branches
uniform_tree(tr)                      # every branch counts equally

Is there more PD than expected?

Raw PD mostly tracks richness. To see more, compare each site to a null expectation.

Null models test hypotheses

Different nulls, different questions

Matrix algorithms run on binary data (fun = "nullcat") or quantitative data (fun = "quantize"); fun = "nullmodel" adds {vegan}’s options.

Randomizing quantitative data

Each site keeps its values; each species keeps its mix of strata. quantize() from {nullcat}

How many iterations?

iters <- ps_suggest_n_iter(
  ps, fun = "quantize",
  method = "curvecat", n_iter = 3e5)

Running a randomization

set.seed(1)
rand <- ps_rand(
  ps, metric = c("PD", "PE", "CE", "RPE"),
  fun = "quantize", method = "curvecat",
  n_rand = 100, n_iter = 5e4)  # ~30 s

Neo- and paleo-endemism (CANAPE)

canape <- ps_canape(rand)

What does phylogenetic turnover measure?

Phylogenetic dissimilarity is the share of branch length not shared by two sites. Species turnover treats every pair of different species as equally different; phylogenetic turnover is low when sites hold close relatives, and high when they hold distant lineages.

Dissimilarity from a focal site

ps_dissim_focal() maps how different every other site is from a focal location. Here, sites in the Sierra Nevada and the Mojave Desert.

focal <- sf::st_as_sf(data.frame(lon = c(-119.6, -116), lat = c(37.8, 34.6)),
                      coords = c("lon", "lat"), crs = 4326)
d_focal <- ps_dissim_focal(ps, focal)

Many kinds of dissimilarity

A wide range of dissimilarity measures are available, including any {vegan} index and the turnover and nestedness components of Sorensen (Baselga & Orme 2012).

sierra <- focal[1, ]
ps_dissim_focal(ps, sierra, method = "sorensen")             # default
ps_dissim_focal(ps, sierra, method = "sorensen_turnover")    # lineage replacement
ps_dissim_focal(ps, sierra, method = "sorensen_nestedness")  # richness differences
ps_dissim_focal(ps, sierra, endemism = TRUE)    # weight range-restricted branches
ps_dissim_focal(ps, sierra, normalize = TRUE)   # proportional composition
ps_dissim_focal(ps, sierra, tips_only = TRUE)   # species only, no tree

An asymmetric index

method = "focal_difference" is the share of the focal site’s branch length missing from each other site. Unlike Sorensen, it ignores lineages the focal site lacks, so it shows where a site’s own lineages occur. For the lineage-poor desert site, they occur almost everywhere.

mojave <- focal[2, ]
d_mojave <- c(ps_dissim_focal(ps, mojave),
              ps_dissim_focal(ps, mojave, method = "focal_difference"))

All pairs of sites

ps <- ps_add_dissim(ps)  # Sorensen

Using the dissimilarity matrix

geo <- ps_geodist(ps)
mt <- vegan::mantel(ps$dissim, geo)
round(c(mt$statistic, mt$signif), 3)  # r, p

[1] 0.406 0.001

Ordination

ps_ordinate() reduces the dissimilarity matrix to a few axes that capture the main gradients in lineage composition (classical or nonmetric multidimensional scaling, or PCA).

ord <- ps_ordinate(ps, method = "cmds", k = 3)

Mapping turnover as color

rgb <- ps_rgb(ps, method = "nmds",
              trans = rank)

Phylogenetic regions

regions <- ps_regions(ps, k = 4,
                      method = "ward.D2")

ps_regions_eval(ps, k = 1:15,
                method = "ward.D2")

What is conservation prioritization?

Choosing where to add new protection, given limited resources.

• Complementarity: a site’s value depends on what’s already protected and what other chosen sites already cover
• Efficiency: maximize biodiversity gained per unit cost
• Targets or benefits: aim to protect a set share of each feature’s range, or maximize total benefit with diminishing returns

In phylospatial, the “features” are branches of the tree, so priorities favor sites that add evolutionary history not yet represented in protected areas.

Prioritization approaches

phylospatial supports two complementary approaches:

• ps_prioritize() ranks every site by its added conservation value, building the ranking one site at a time. The result is nested: any top-k sites form a sensible plan, which is useful when budgets are uncertain or protection will happen incrementally.
• ps_prioritizr() sets up a problem for the {prioritizr} package, with a range-protection target or a budget for every branch, and an integer programming solver finds the optimal set of sites for that one scenario.

Prioritization inputs

Existing protection (init, from 0 to 1) and the relative cost of protecting each site, both on the moss grid (ext is the package data folder from earlier).

init <- rast(file.path(ext,
                       "moss_protection.tif"))

pop  <- rast(file.path(ext,
                       "moss_popdens.tif"))
cost <- 1 + log1p(pop)

Diminishing returns

plot_lambda()

A priority ranking

priority <- ps_prioritize(
  ps, lambda = 1, init = init, cost = cost)

How much does each step buy?

perf <- ps_performance(ps, priority,
                       target = 0.8)
plot(perf, yvar = "cov80")

Optimal solutions with {prioritizr}

library(prioritizr)
prob <- ps_prioritizr(
  ps, init = init, cost = cost,
  objective = "min_set", target = 0.3)
solution <- solve(prob)

Flexible inputs. One object. Many analyses.

Thank you!