AlleleShift: an R package to predict and visualize population-level changes in allele frequencies in response to climate change

Data import

Genetic response data (including a matrix with populations as rows and allele counts as columns) for the calibration of the AlleleShift::count.model and prediction via AlleleShift::count.pred is required to be in the adegenet::genpop format. Individual-based data that are in the adegenet::genind format can be converted into the genpop format via the adegenet::genind2genpop function. The adegenet and poppr packages provide various methods of importing data from other software application formats into the genind format, such as adegenet::import2genind and poppr::read.genalex. Environmental data of populations, used as explanatory variables in redundancy analysis (RDA), is expected to be provided as a data.frame with the same sequence of populations as the genetic response data (this is a general requirement for community ecology methods in the vegan and BiodiversityR packages; a check is available via BiodiversityR::check.datasets). Whereas environmental data typically is baseline and changed (bio)climatic data such as is available from WorldClim (Fick & Hijmans, 2017), ENVIREM (Title & Bemmels, 2018), CHELSA (Karger, Schmatz & Dettling, 2020) or PaleoClim (Brown et al., 2018), it is also possible to expand input to other data available for species distribution modelling (a recent overview of available data sets is provided by Booth (2018).

Data analysis

Prior to calibrating the models that predict allele frequencies from bioclimatic explanatory data, it is recommended to reduce the explanatory variables to a subset where the Variance Inflation Factor (VIF) is below a predefined threshold for each variable. Such methods are also recommended for regression analysis (Fox & Monette, 1992) and species distribution modelling (Kindt, 2018). With this approach, it is easy to select the same variables from future data sets and comparison with other studies may also become easier. VIF analysis, and an additional feature of forcing the algorithm to keep preselected variables within the final subset, is available via AlleleShift::VIF.subset (Table 1). This is a function that uses BiodiversityR::ensemble.VIF.dataframe internally after a first step of removing all explanatory variables that have correlations larger than the VIF threshold with the preselected variables. There is an option to generate a correlation matrix for the final subset of variables (Fig. 1). I also recommend conducting the VIF analysis for the genetic response data (see discussion and Fig. 1). I further advise to remove any individuals with partially missing genetic or (bio)climatic data prior to the analysis. As a sensitivity analysis, results could possibly be compared with models calibrated with data where missing genotypes were replaced by genotypes predicted by the impute function of the LEA package (Gain & François, 2021).

Prior to model calibration, I suggest checking (according to criteria described below) for the shifts of populations in environmental space between the baseline and changed climates. This can be achieved via function AlleleShift::population.shift, which draws arrows between each population in the baseline and changed climate. There are alternatives of using principal components analysis (PCA) or redundancy analysis (RDA) to generate the ordination diagrams (Fig. 2). What is important to check in the ordination graphs is whether populations shift in a similar fashion, as that will facilitate the interpretation of predicted shifts in allele frequencies. If some populations would show a different trend in the ordination graph, their allele frequencies would also be expected to change in a different way.

With the selected genetic and explanatory data, model calibration can be done. In a first step, a RDA model is fitted (AlleleShift::count.model) that can predict counts of alleles in baseline or changed climates (AlleleShift::count.pred). The user has the option also to obtain results for the canonical correspondence analysis procedure of Blumstein et al. (2020) with the count model via argument cca.model=TRUE. In the second step, the predicted allele counts for the baseline climate serve as explanatory variables for the calibration of a generalized linear model (GAM via mgcv::gam; Wood, 2004) with the baseline allele frequencies as response and a binomial family function (AlleleShift::freq.model). This procedure ensures that predictions are within the realistic interval for frequencies between 0 to 1. Function AlleleShift::freq.pred allows the prediction of allele frequencies for baseline and changed climates. A simulation study (Article S2, using a subset of a simulation study available from Frichot & François (2015)) and the results below show that the predicted frequencies approximate the expected frequencies well.

The output of the two steps (RDA followed by GAM) is a data.frame as shown in Table 2 (to fit printing space available here, the data is shown in a transposed format where rows and columns were swapped). All figures shown in this manuscript were obtained with the example data sets of AlleleShift::Poptri.genind (individually-based allele counts), AlleleShift::Poptri.baseline.env (climatic descriptors of the populations in the baseline climate), AlleleShift::Poptri.future.env (climatic descriptors of the populations in the future climate) and AlleleShift::Poptri.loc (geographical coordinates of the populations). These data sets were converted from the example data sets provided by Blumstein et al. (2020) for Populus trichocarpa. It can be seen for population Nisqually in our case study that negative allele counts and frequencies are predicted for one of the minor alleles in the RDA prediction step, but that the GAM step predicts the biologically acceptable frequency of 0.027. Function AlleleShift::freq.ggplot (Fig. 3) enables a visual inspection of the power of the models to predict allele frequencies for the calibration data.

Visualizations

AlleleShift generates various types of ggplot2 (Wickham, 2016; version 3.3.2) graphics from the output of AlleleShift::freq.pred. These include dot (Fig. 4), pie or donut (Fig. 5), moon (Fig. 6) and waffle (Fig. 7) graphics, and smoothed regression surfaces (Fig. 8). As an intermediate step to generate various of these graphics, helper functions such as waffle.baker for shift.waffle.ggplot or moon.waxer for shift.moon.ggplot prepare data for the main graphing function (Table 1). With default settings, visualizations depict changes in allele frequencies for each allele in different panels, internally via ggplot2::facet_grid. Setting argument mean.change to TRUE, the graphics depict median or mean changes in allele frequencies.

Function shift.surf.ggplot plots populations in geographical space via their geographical coordinates (longitude and latitude in Fig. 8). The function fits a smoothed regression surface for allele frequencies via vegan::ordisurf. This output is then further processed internally within the function via BiodiversityR::ordisurfgrid.long (an overview of generating ggplot2 ordination diagrams via vegan and BiodiversityR is given in Kindt (2020b); these guidelines can be used to customize ordination graphs as shown in Figs. 2 and 8). Various options of fitting smoothed regression surfaces are available by providing additional arguments to shift.surf.ggplot, such as the various mgcv::smooth.terms options of thin plates, Duchon splines, cubic regression splines, P-splines, Markov Random Fields, Gaussian process smooths, soap film smooths, splines on the sphere and random effects.

As graphics are generated with ggplot2, it is relatively easy to generate animated versions of visualizations with gganimate (Pedersen & Robinson, 2020; version 1.0.7). Scripts for generating animated versions for dot graphics, pie graphics and smoothed surfaces are provided with the documentation of the respective functions. Video S3–S5 are for a time series that interpolates bioclimatic data between baseline and future climate in steps of five years; these also extrapolate data far into the future (up to a million years, see the Results and Discussion for my reason to do this). Other than time series, animated visualizations could also be generated for different global circulation models (GCMs) or scenarios such as various shared socio-economic pathways (SSPs) developed in the context of the sixth assessment report of the Intergovernmental Panel on Climate Change.

Results and Discussion

AlleleShift predicts shifts in allele frequencies via RDA and GAM, an alternative pathway that maintains Euclidean distances among populations and individuals (Article S1). It also avoids making negative frequency predictions as what may occur with the protocol recently devised by Blumstein et al. (2020). My methodology however faces the same limitations of data requirements as discussed by Blumstein et al. (2020) for their protocol in terms of the initial identification of responsive markers. Key assumptions of the Blumstein et al. (2020) protocol apply also to my approach, including that allelic effects are independent and additive with no epistasis or dominance (see also the discussion on epistasis, structural genomic variation and epigenetics by Stange, Barrett & Hendry (2020).

I recommend reducing both the explanatory variables and the response variables to subsets of data with a maximum Variance Inflation Factor (VIF) of 20 or less (Ter Braak & Smilauer, 2002), for alleles to obtain better estimates of changes in their frequencies. For the allele counts as response variables, if VIF would be larger than 20, I would use function VIF.subset in an iterative procedure whereby earlier selected subsets of alleles are excluded, to generate different subsets of alleles (but keeping variables for both A and B allele counts in each final subset). For example, with 20 alleles X01A, X01B, X02A, X02B, …, X10A, X10B, the two subsets could be {X01A, X01B, …, X05A, X05B} and {X06A, X06B, …, X10A, X10B}. With the various subsets, shifts in allele frequencies can then be predicted, and finally predictions with all subsets can be combined.

When predictions are made by AlleleShift into the future, and especially into novel climatic conditions, it is warranted to consider the transferability of the calibrated models and ideally to provide ‘transferability metrics’ that quantify prediction uncertainty (Yates et al., 2018). The function of environmental.novel identifies populations that have novel climatic conditions (values outside the baseline range) for at least one of the explanatory variables for the changed climate. This function also calculates the probabilities of encountering the novel conditions from the mean and standard deviation of the baseline ranges and then returns the lowest of these probabilities as an indicator of the degree of novelty of the environmental conditions.

For the animated graphics, as an extreme example of the caveats of the methodology, I made projections one million years into the future. These projections followed earlier trends of the 21st century and resulted in allele frequencies becoming fixed at 0 or 1, which is a biologically possible scenario. At the same time, however, these simulations clearly illustrated that the methodology through RDA is correlative as in correlative/phenomenological species distribution models. In fact, the linear extrapolation of climate variables resulted in an environmental data set where the mean annual temperature was above 30,000 degrees, which obviously is a biologically irrelevant scenario. Thus, the method presented here should be used cautiously in novel climates especially as predictions will be available (the model will not crash or give an error warning), and should thus be used even more cautiously where differences between baseline and future climates are large.

As an approach to obtain a better handle on the transferability of AlleleShift predictions into future climates, I recommend to also estimate the habitat suitability via species distribution models (SDMs). Where habitat suitability is below a certain threshold, for instance a threshold where the sum of sensitivity and specificity is maximized (Liu, White & Newell, 2013), results of predictions of AlleleShift could be annotated of having lower transferability. Similarly, where evaluation strips as proposed by Elith et al. (2005) suggest that it is less likely that the habitat is suitable within the changed range of environmental conditions, additional information could be provided on a lower transferability. Well-documented methods of utilizing SDMs to predict shifts in species habitat suitability are available in the literature, including recent examples that use the ensemble suitability modelling framework available in BiodiversityR (e.g., de Sousa et al., 2019; Fremout et al., 2020; Kindt, 2018; Ranjitkar et al., 2014). For organisms such as trees, correlative SDMs remain the best available method of predicting future species suitability, whereas the limitations of these methods may not be as great as has been suggested (Booth, 2016; Booth, 2018). What is also attractive about SDMs is that a wider set of presence observations are likely to be available than those populations that have been studied genetically. Presence data are available from open-source databases such as GBIF or the Botanical Information and Ecology Network (Enquist et al., 2016). Further to the collation of a larger set of presence observations, application of SDMs should be straightforward using the same (bio)climatic data sets as applied in AlleleShift. The approach of combining the results of AlleleShift with SDM is somewhat similar to the method applied by Aguirre-Liguori et al. (2019) to develop species distribution models for alleles. In my proposal, however, the predictions of allele frequencies and SDM are done independently, and ideally with an expanded point presence data set for SDM.

A straightforward and practical expansion of the methodology I have proposed is to tree seed sourcing programmes (Broadhurst et al., 2008), possibly for developing schemes of human-assisted geneflow sensu Aitken & Whitlock (2013). This is important for ensuring the matching of planting materials to the conditions of planting sites (Cernansky, 2018; Roshetko et al., 2018; Kettle et al., 2020). For specific planting sites and planting times (considering the perennial nature of trees, climate change during the production cycle should be considered) of interest, the prediction methods can readily provide the predicted allele frequencies needed for adaptation. Theoretically, based on the similarity between predicted allele frequencies and those of available source populations, the best matching source can then be selected. Similar approaches to devise transfer and conservation schemes in the face of climate change can be employed for other organisms than plants (Fitzpatrick & Edelsparre, 2018; Rochat, Selmoni & Joost, 2021).