Genotype–environment association in Swiss stone pine

0. Technical details

• Similar to last time, you likely need to install some packages as outlined below.
• Navigate to an easy to find directory in R, potentially the one you make last time, with setwd("path/to/your/directory")
• Download the data from drive

1. The dataset and the question

Pinus cembra grows at the Alpine timberline. Alpine plants have been intensively studied for strong selective pressures such as frost tolerance, high UV exposure, and the ecological differences over relatively small distances. Dauphin et al. (2020) sequenced the exome of 480 trees from 24 populations across the Swiss Alps and asked which SNPs are associated with environment more than population structure. If there are signals of local adaptation among some populations, that could be important for timber management.

Here, we have an allele frequency matrix (populations × sites) and an environmental matrix (populations × variables). We will use this to find SNPs that covary with environment (in a multideimensional sense). The difficulty is that populations close in space can give associations with environment due to shared ancestry rather than local adaptation. We will use the redundancy analysis (RDA) and related Genotype-Environment Association (GEA) methods for separating those two components.

Two caveats about this data:

• The genotypes are population allele frequencies, not individual genotypes. This is like estimating \(p\) from the population.
• There are 24 populations. That is your sample size since we pool by population, not 480. We can discuss why you might do this for research.

2. Reading the data into R

library(vegan)      # rda(), varpart(), ordiR2step(), anova.cca()
library(robust)     # covRob() -- robust covariance for the genome scan
library(lfmm)       # lfmm_ridge(), lfmm_test() -- the paper's own method
library(ggplot2)
library(dplyr)
library(corrplot)
library(ggrepel)
library(patchwork)

# use 9999 for anything you would publish.
NPERM <- 999
set.seed(973479)

Prepare the allele frequency data for the ordination analyses:

snp <- read.csv("data/Pinus_cembra_All_17061_SNPs_MEC-19-1345.csv")
env <- read.csv("data/Pinus_cembra_All_Env_variables_MEC-19-1345.csv")

# Columns 1-2 are CHROM and POS; the other 24 are populations. Transpose so that
# populations are rows, which is what the downstream ordination function in R expects.
Y <- t(as.matrix(snp[, -(1:2)]))

# CHROM repeats across SNPs on the same contig, so CHROM_POS is the unique id
colnames(Y) <- paste(snp$CHROM, snp$POS, sep = "_")

dim(Y)
Y[1:4, 1:5]

Before anything else, we need to check that the 24 populations are in the same order for both the genetic and environmental data. The analyses will not tell you this, it is our responsibility to ensure the data is organized correctly.

# If this ever failed we would be pairing one population's genetics with
# another's climate, and every result below would be silently wrong.
stopifnot(identical(rownames(Y), env$Sample_ID))
cat("population IDs match:", nrow(Y), "populations\n")
cat("allele frequency range:", range(Y), " missing values:", sum(is.na(Y)), "\n")

There is no missing data at all, a benefit of working with populations rather than individuals. Now, we filter away sites with very low or very high minor allele frequency. Those are near fixation and likely not very informative about the task at hand.

# A SNP that is nearly fixed everywhere carries almost no information about
# environment, and its frequencies are mostly sampling noise. Drop both tails.
freq_mean <- colMeans(Y)
Yf <- Y[, freq_mean > 0.05 & freq_mean < 0.95]
cat(ncol(Y), "SNPs ->", ncol(Yf), "after MAF filtering\n")

You should be left with 12,506 SNPs.

3. The environmental variables

There are 34 predictors: 19 standard bioclim variables (mostly about temperature and precipitation) and 15 topographic ones derived from a digital elevation model (DEM). These include things like slope, and drainage. The topographic set is the point of this study system. In the mountains, two populations 5 km apart can share a climate and differ in topographic features.

bio_vars  <- paste0("bio", 1:19)                              # climate
topo_vars <- grep("^t[0-9]{2}_", names(env), value = TRUE)    # terrain

# The environmental variables need to be standardised (centered and scaled) so that each has a mean of 0 and a standard deviation of 1. This prevents variables with larger scales (e.g. altitude) from dominating the analysis.
E <- scale(env[, c(bio_vars, topo_vars)])
cat(length(bio_vars), "bioclim +", length(topo_vars), "topographic =", ncol(E), "predictors\n")

We can now check the correlations among the environmental variables:

corrplot(cor(E[, bio_vars]), type = "upper", tl.cex = 0.7, tl.col = "black")

Pairwise correlations among the 19 bioclim variables. Blue is positive, red negative; the darker the square the stronger the correlation.

This is a common problem. Multiple Bioclim variables are summaries of the same monthly temperature and precipitation data. It is nosurprise that there is some autocorrelation.

cm <- cor(E)
diag(cm) <- 0                                  # ignore each variable's correlation with itself
hi <- which(abs(cm) > 0.95, arr.ind = TRUE)    # row/column indices of the strong pairs
hi <- hi[hi[, 1] < hi[, 2], , drop = FALSE]    # keep each pair once, not twice
data.frame(var1 = rownames(cm)[hi[, 1]], var2 = colnames(cm)[hi[, 2]],
           r = round(cm[hi], 3))

Nine pairs correlate above 0.95, topping out at 0.995 for bio14 (precipitation of the driest month) against bio17 (precipitation of the driest quarter). Using both in a model would not be helpful. A common approach is to prune correlated variables iteratively until nothing exceeds a threshold.

prune_correlated <- function(E, threshold) {
  keep <- colnames(E)                   # start with every variable
  repeat {
    cm <- cor(E[, keep])
    diag(cm) <- 0                       # a variable always correlates 1 with itself
    if (max(abs(cm)) < threshold) break # nothing left above the threshold: done

    # locate the single most correlated remaining pair
    w <- which(abs(cm) == max(abs(cm)), arr.ind = TRUE)[1, ]

    # Of that pair, discard whichever is more redundant with everything else,
    # so we keep the more independent member.
    mean_cor_1 <- mean(abs(cm[w[1], ]))
    mean_cor_2 <- mean(abs(cm[w[2], ]))
    drop_it <- if (mean_cor_1 > mean_cor_2) keep[w[1]] else keep[w[2]]

    keep <- setdiff(keep, drop_it)
  }
  keep
}
for (thr in c(0.9, 0.8, 0.7)) {
  cat(sprintf("|r| < %.1f  ->  %2d predictors\n", thr, length(prune_correlated(E, thr))))
}
keep_vars <- prune_correlated(E, 0.7)
keep_vars

The 0.7 threshold is commonly cited and leaves 14 predictors. This is a judgement call, since a different threshold would affect the predictors in your model.

4. Population structure

Before we constrain anything, it is worth being clear about what these methods are doing, because RDA sounds more exotic than it is.

Yf is a table of 24 populations by 12,506 SNPs. You cannot look at 12,506 dimensions, so you need a way to draw it in some dimensional space we can interpret. That means two or three at most. Ordination finds these axes for us, using weighted combinations of the original columns that maximize the amount of variation explained.

PCA asks: what are the dominant axes of variation in this table? It only considers the genetic data.

RDA asks a more specific question: what are the dominant axes of the variation that climate can predict? There are two steps here. First, it regresses every SNP on the climate variables and keeps the fitted values — the part of each SNP’s frequencies that climate accounts for. Then it runs a PCA on those fitted values. So an RDA really is a PCA, just of a filtered version of the data.

That is where “constrained” comes from. In a PCA the axes can point anywhere. In an RDA they are required to be linear combinations of the predictors, so there can never be more constrained axes than predictors — with four climate variables we get four RDA axes, no matter how many SNPs we have.

Run with no predictors at all, rda() does an ordinary PCA, and that is how we summarise population structure here.

pca <- rda(Yf, scale = FALSE)
eig <- pca$CA$eig
round(100 * eig[1:6] / sum(eig), 1)

# "sites" = the rows of the response matrix, i.e. our 24 populations
sites <- as.data.frame(scores(pca, choices = 1:2, display = "sites", scaling = 0))
sites$pop    <- rownames(Yf)
sites$region <- ifelse(grepl("^CH_", sites$pop), "CH", "HJ")

ggplot(sites, aes(PC1, PC2, colour = region)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey80") +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey80") +
  geom_point(size = 2.6) +
  geom_text_repel(aes(label = pop), size = 2.6, max.overlaps = 20, show.legend = FALSE) +
  labs(x = sprintf("PC1 (%.1f%%)", 100 * eig[1] / sum(eig)),
       y = sprintf("PC2 (%.1f%%)", 100 * eig[2] / sum(eig)), colour = NULL) +
  theme_bw(base_size = 11)

Genetic PCA of the 24 populations.

PC1 carries 15.5% and the first three together 29.3%. These are wind-pollinated conifers, so similar to our big bluestem, we don’t really expect large variances explained by a single axis. We keep the first three PCs as a decent approximation of the population structure.

Now build the predictor table every model from here on will use: the 34 standardised environmental variables, plus the first three genetic PCs as a summary of shared ancestry.

vars <- data.frame(E)
vars <- cbind(vars, as.data.frame(scores(pca, choices = 1:3, display = "sites", scaling = 0)))
names(vars)[(ncol(vars) - 2):ncol(vars)]

To see what “constrained” buys us, ordinate with a few environmental variables and compare to the unconstrained PCA. We will pick the climate variables next but are peeking ahead to make a point.

# the constrained ordination, using the variables section 5 will select
rda_peek <- rda(Yf ~ bio4 + bio18 + bio8 + bio3, data = vars)

# pull the population coordinates out of each ordination separately, so that
# both tables have the same three columns before we stack them
ord_pca <- as.data.frame(scores(pca, choices = 1:2, display = "sites", scaling = 0))
colnames(ord_pca) <- c("axis1", "axis2")
ord_pca$pop   <- rownames(Yf)
ord_pca$panel <- "PCA (unconstrained)"

ord_rda <- as.data.frame(scores(rda_peek, choices = 1:2, display = "sites", scaling = 0))
colnames(ord_rda) <- c("axis1", "axis2")
ord_rda$pop   <- rownames(Yf)
ord_rda$panel <- "RDA (constrained by climate)"

ord <- rbind(ord_pca, ord_rda)
ord$region <- ifelse(grepl("^CH_", ord$pop), "CH", "HJ")

ggplot(ord, aes(axis1, axis2, colour = region)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey85") +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey85") +
  geom_point(size = 2.2) +
  geom_text_repel(aes(label = pop), size = 2.2, max.overlaps = 30, show.legend = FALSE) +
  facet_wrap(~ panel, scales = "free") +
  labs(x = "axis 1", y = "axis 2", colour = NULL) +
  theme_bw(base_size = 11)

The same 24 populations ordinated without constraint (left) and constrained by four climate variables (right). Very little change.

# how similar are the two first axes? (1 would mean identical ordering)
cor(ord_pca$axis1, ord_rda$axis1)

The two panels are very similar, and the first axes correlate at 0.994.

This shows the frequent confounding of structure and other factors we are interested in. We can fix this though!

Q1: If climate were shaping genetic variation independently of shared ancestry, what would you expect this PCA-vs-RDA comparison to look like instead? How does associations between genetic variation and environment arise without local adaptation?

5. Choosing predictors

We need a model with predictors (e.g. envioronment and topography). More predictors are better, right?… Try to use ordiR2step on all 34 variables:

full_all <- rda(as.formula(paste("Yf ~", paste(c(bio_vars, topo_vars), collapse = " + "))),
                data = vars)
RsquareAdj(full_all)$adj.r.squared

\(R^2\) = NA. With 24 populations and 34 predictors, the model has more parameters than observations, eleven terms get aliased, and adjusted R² is undefined. ordiR2step is nice that it stops rather than giving a bad result, but not all programs handle this gracefully.

Now the same procedure on the pruned set after the iterative elimination:

full <- rda(as.formula(paste("Yf ~", paste(keep_vars, collapse = " + "))), data = vars)
null <- rda(Yf ~ 1, data = vars)
# R2adj = adjusted R-squared: the share of genetic variation the predictors explain,
# penalised for how many predictors were used. With n = 24 that penalty is severe.
cat("full model R2adj:", round(RsquareAdj(full)$adj.r.squared, 4), "\n")

step <- ordiR2step(null, full, Pin = 0.05, R2permutations = 500,
                   permutations = NPERM, trace = FALSE)
step$anova

Four variables are retained (at least on my last run): bio4 (temperature seasonality), bio18 (precipitation of thewarmest quarter), bio8 (mean temperature of the wettest quarter) and bio3 (isothermality). All temperature-regime and summer-moisture variables, plausible for a timberline conifer. Notably, no topographic variable survives, despite the study being designed around topography. This frequently happens when people aim to do interesting research - temperature and precipitation still wins out.

A stable selection procedure should recover the same predictors every time, given the same data and the same set.seed(). That is not always guaranteed in practice: ordiR2step chooses variables using permutation tests, and with n = 24 the exact permutation draws might affect the predictors here. I have not thoroughly evaluated here.

sel_vars <- setdiff(all.vars(formula(step))[-1], "Yf")
sel_vars
# formulas below spell these four out literally; edit them if your run selects different variables
round(vif.cca(rda(Yf ~ bio4 + bio18 + bio8 + bio3, data = vars)), 2)

Variance inflation factors (VIFs) below 5 are comfortable; these are all well under.

Q2: No topographic variable made it into the selected model, even though the arguments for including them are compelling. What is one statistical reason that could happen, and what is a separate biological reason? Can you think of a simple analysis that would help you tell those two explanations apart?

6. RDA, partial RDA, and where the signal actually lives

Now the two models that matter. The first regresses allele frequencies on climate. The second does the same while conditioning on the genetic PCs.

Condition() is worth understanding. It fits the conditioning variables first, takes the residuals the difference between the observed and fitted values. This conceptually translates to what is left of the SNP frequencies once population structure has been subtracted, and hands only the remaining variation in allele frequencies to the climate model. Climate is then credited with what structure could not already explain.

This splits the total variation three ways, and the next chunk prints all three: conditioned (the PCs’ share), constrained (climate’s share of what was left) and residual (neither).

# climate only
f_rda  <- Yf ~ bio4 + bio18 + bio8 + bio3

# climate, after population structure has been given its share first
f_prda <- Yf ~ bio4 + bio18 + bio8 + bio3 + Condition(PC1 + PC2)

mod_rda  <- rda(f_rda,  data = vars)
mod_prda <- rda(f_prda, data = vars)

# tot.chi is total variation; CCA is the constrained part, pCCA the conditioned part

cat("RDA   R2adj:", round(RsquareAdj(mod_rda)$adj.r.squared, 4),
    "  constrained:", sprintf("%.1f%%", 100 * mod_rda$CCA$tot.chi / mod_rda$tot.chi), "\n")
cat("pRDA  R2adj:", round(RsquareAdj(mod_prda)$adj.r.squared, 4),
    "  constrained:", sprintf("%.1f%%", 100 * mod_prda$CCA$tot.chi / mod_prda$tot.chi),
    "  conditioned:", sprintf("%.1f%%", 100 * mod_prda$pCCA$tot.chi / mod_prda$tot.chi), "\n")

Climate explains R²adj ≈ 0.150 on its own. Condition on two genetic PCs and that collapses to ≈ 0.038. Three-quarters of the apparent environmental signal was equally well explained by population structure.

Both models are still significant, so the environmental effect is real, not nothing. With NPERM <- 999 the smallest p-value obtainable is 1/1000 = 0.001, so a reported 0.001 is still meaningful.

anova.cca(mod_prda, permutations = NPERM)

Variance partitioning splits total variation into what only climate explains, what only structure explains, what both explain, and what neither does.

# Splits total variation four ways: climate only, structure only, the part both
# explain equally well, and the part neither explains.
vp <- varpart(Yf, ~ bio4 + bio18 + bio8 + bio3, ~ PC1 + PC2, data = vars)
vp

plot(vp, digits = 3, Xnames = c("Climate", "Structure"), bg = c("#F9A242", "#6B4596"))

Variance partitioning. [a] climate alone, [b] population structure alone, [c] the fraction they share, [d] residual.

Q3: The shared fraction here (0.112) is larger than either the climate-only or structure-only fraction. What does that ordering suggest about how confidently you can attribute genetic variation to climate versus shared ancestry in this dataset? Can you think of a sampling design that would shrink this confounded fraction?

7. The genome scan

Fitting an RDA tells you climate matters overall. It does not tell you which SNPs carry the signal. For that we need the loadings.

A loading is a SNP’s coordinate on a constrained axis. Recall from section 4 that each RDA axis is a particular weighted combination of the climate variables — RDA1 here is mostly a temperature- seasonality contrast. A SNP’s loading on RDA1 says how strongly its allele frequencies move along that climate gradient across the 24 populations. Near zero means the SNP’s frequencies are unrelated to it; large positive or negative means they track it closely, in one direction or the other.

If you plot every SNP by its loadings on the first two axes you get a dense cloud centred on the origin. Most SNPs are in that blob: their frequencies vary for reasons that have nothing to do with climate. The ones worth a second look are the ones far from the centre in any direction — and that is a geometric question, which makes the genome scan an outlier-detection problem rather than a significance test on each SNP separately.

Distance from the centre of a cloud is measured with a Mahalanobis distance, which accounts for the axes having different spreads and being correlated. A normal Euclidean distance would over-weight whichever axis happens to be widest. And it has to be a robust Mahalanobis distance: the covariance has to be estimated from the same cloud we are hunting outliers in, so a handful of genuinely extreme SNPs would inflate the estimated spread and make themselves look ordinary. covRob() estimates this appropriately.

One choice remains: K, how many constrained axes count toward the distance. We use 2, meaning a SNP must be unusual with respect to the first two climate axes. In practice, you might want to explore the stability of candidate SNPs across different choices of K.

# Adapted from Capblancq & Forester's rdadapt(). One change: the original uses
# Storey q-values; we use Benjamini-Hochberg, which is base R and is also what
# Dauphin et al. used (BH at FDR 0.01).
rdadapt_bh <- function(mod, K) {

  # Each SNP's coordinates on the first K constrained axes.
  loadings <- mod$CCA$v[, 1:K, drop = FALSE]

  # Distance of every SNP from the centre of the loading cloud. Mahalanobis
  # because the axes have different spreads; robust (covRob) because the spread
  # is estimated from the very cloud we are hunting outliers in, and a few
  # extreme SNPs would otherwise inflate it and hide themselves.
  D <- covRob(scale(loadings), distance = TRUE,
              na.action = na.omit, estim = "pairwiseGK")$dist

  # Rescale so the median distance matches the chi-square expectation. This is
  # the same genomic-inflation idea used by the LFMM analysis of Dauphin et al. 2020
  lambda <- median(D) / qchisq(0.5, df = K)

  # Under the null, the rescaled distance follows chi-square with K df
  p <- pchisq(D / lambda, K, lower.tail = FALSE)

  # Benjamini-Hochberg: control the share of false positives among the SNPs we call
  data.frame(snp     = rownames(mod$CCA$v),
             p.value = p,
             q.value = p.adjust(p, method = "BH"),
             row.names = NULL)
}

scan_prda <- rdadapt_bh(mod_prda, K = 2)
cand_prda <- scan_prda$snp[scan_prda$q.value < 0.05]
cat("candidate SNPs (pRDA, BH q < 0.05):", length(cand_prda),
    sprintf("(%.2f%% of %d)", 100 * length(cand_prda) / ncol(Yf), ncol(Yf)), "\n")
head(scan_prda[order(scan_prda$q.value), ], 5)

man <- data.frame(idx  = seq_len(nrow(scan_prda)),
                  p    = scan_prda$p.value,
                  type = ifelse(scan_prda$q.value < 0.05, "Candidate", "Neutral"))
man <- man[order(man$type == "Candidate"), ]   # candidates drawn on top of the grey

ggplot(man, aes(idx, -log10(p), colour = type)) +
  geom_point(size = 0.9) +
  scale_colour_manual(values = c(Neutral = "grey85", Candidate = "#F9A242")) +
  labs(x = "SNP (input order)", y = expression(-log[10](p)), colour = NULL) +
  theme_bw(base_size = 11)

Genome scan. Each point is one SNP; candidates pass Benjamini-Hochberg at FDR 5%.

The x-axis is input order, not genomic position — this is a de-novo transcriptome assembly with 4,677 contigs, so there is no chromosome to plot along. Worth keeping in mind that runs of adjacent points here don’t reflect linkage.

# vegan calls the columns of the response matrix "species"; here they are SNPs
loc <- as.data.frame(scores(mod_prda, choices = 1:2, display = "species", scaling = "none"))
loc$type <- ifelse(rownames(loc) %in% cand_prda, "Candidate", "Neutral")
loc <- loc[order(loc$type == "Candidate"), ]   # draw candidates last, on top

# "bp" = biplot scores: where each climate variable points in this space
arrows_df <- as.data.frame(scores(mod_prda, choices = 1:2, display = "bp"))

ggplot() +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "grey80", linewidth = 0.6) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey80", linewidth = 0.6) +
  geom_point(data = loc, aes(RDA1 * 20, RDA2 * 20, colour = type), size = 1.2) +
  scale_colour_manual(values = c(Neutral = "grey85", Candidate = "#F9A242")) +
  geom_segment(data = arrows_df, aes(x = 0, y = 0, xend = RDA1, yend = RDA2),
               arrow = arrow(length = unit(0.02, "npc")), colour = "black", linewidth = 0.4) +
  geom_text(data = arrows_df, aes(1.12 * RDA1, 1.12 * RDA2, label = rownames(arrows_df)),
            size = 3) +
  labs(x = "RDA 1", y = "RDA 2", colour = NULL) +
  theme_bw(base_size = 11)

SNP loadings in RDA space (scaled up 20x for visibility) with climate variables as arrows. Candidates sit far from the origin in the direction of the variables they track.

Q4: A SNP earns “candidate” status here just by sitting far from the centre of the loading cloud, estimated from only 24 populations. What would make you more confident that a given candidate reflects real selection rather than sampling noise? What would make you less confident?

8. Takeaways

1. Your sample size is the number of populations. 480 trees and 17,061 SNPs sound like a lot; n = 24 is what the model is trying to fit, and it is why the 34-predictor model was not possible.
2. Collinear predictors are the normal and we deal with them. Bioclim variables are derived from each other. Nine pairs here exceed |r| = 0.95. Prune before you model.
3. Most of the climate signal was confounded with ancestry. The shared fraction (0.112) is about three times either pure fraction. Conditioning cut climate’s R²adj from 0.150 to 0.038.
4. A residual signal survived. The partial RDA is still significant (p = 0.001 — the smallest value 999 permutations can return). Confounding makes evaluating environmental effects on genetic variation tricky, it does not make the environmental effect disappear.
5. Candidate lists are soft. They move with K, with conditioning, and with method. Treat them as ranked hypotheses for follow-up, and give most weight to SNPs that survive more than one approach. Experimental validation afterwards is ideal but often time- and budget-limited.