Heredity (2009) 102, 330–341

& 2009 Macmillan Publishers Limited All rights reserved 0018-067X/09 $32.00
www.nature.com/hdy
REVIEW

Genetic markers in the playground of

multivariate analysis
T Jombart, D Pontier and A-B Dufour
Université de Lyon, F-69000, Lyon, Université Lyon 1, CNRS, UMR5558, Laboratoire de Biométrie et Biologie Evolutive, Villeurbanne,
France

Multivariate analyses such as principal component analysis review, we provide a critical analysis of the application of
were among the first statistical methods employed to multivariate methods to genetic markers, using a general
extract information from genetic markers. From their early framework that unifies all these methods for the sake of
applications to current innovations, these approaches clarity. First, we focus on some common mistakes in these
have proven to be efficient for the analysis of the genetic applications and ways to avoid these pitfalls. We then detail
variability in various contexts such as human genetics, the most critical particularities of allele frequencies that
conservation and adaptation studies. However, because demand adaptations of multivariate methods, and we
multivariate analysis is a wide and diversified area of propose solutions to the subsequent problems. Finally, we
statistics, choosing a method appropriate to both the data tackle several questions of interest in which multivariate
and to the question being asked can be difficult. Moreover, analysis has a great role to play, such as the study of the
some particularities of genetic markers need to be taken into typological coherence of different genetic markers, or the
account when using multivariate methods. As a conse- investigation of spatial genetic patterns.
quence, multivariate analyses are often used as black boxes, Heredity (2009) 102, 330–341; doi:10.1038/hdy.2008.130;
which results in frequent mistakes in the literature. In this published online 21 January 2009

Keywords: multivariate analysis; ordination; genetic markers; principal component analysis; statistics; methods

Introduction adaptations (Johnson et al., 1969; Mulley et al., 1979;

Barker et al., 1986).
Statistical methods have long become an essential Multivariate analysis has several advantages over
component of the toolbox of population geneticists other classical approaches used in population genetics,
(Fisher, 1952). Developments in statistical theories and like the Bayesian clustering implemented in the software
the continual increases in cheap computing power STRUCTURE (Pritchard et al., 2000; Falush et al., 2003).
provide numerous tools for genetic marker analysis, First, multivariate methods are exploratory, that is, they
allowing geneticists to address new and challenging do not require strong assumptions about an underlying
questions. Multivariate analyses (also called ordinations in genetic model, such as the Hardy–Weinberg equilibrium
reduced space) such as principal component analysis or the absence of linkage disequilibrium. Although
(Pearson, 1901) have been shown to be efficient in clustering approaches suppose that genotypes are
extracting information from genetic markers (Cavalli- structured in discrete populations, ordinations in re-
Sforza, 1966; Johnson et al., 1969; Smouse et al., 1982) duced space simply aim at summarizing the genetic
because of their ability to summarize multivariate variability, and can therefore reveal any kind of genetic
genetic information into a few synthetic variables. From structuring including clines (for example, Jombart et al.,
these early applications to current innovative develop- 2008). Multivariate methods are not computer-intensive,
ments (Patterson et al., 2006; Pavoine and Bailly, 2007; and can be applied to huge datasets (such as ‘hundreds of
Jombart et al., 2008), these methods have proven to be thousands of markers and thousands of samples’ in Patterson
useful in various fields, such as human genetics et al. (2006)), for which Bayesian clustering would be
(Menozzi et al., 1978; Bertranpetit and Cavalli-Sforza, impractical. Moreover, multivariate analysis can address
1991; Cavalli-Sforza et al., 1993), conservation (Moazami- complex questions such as identification of adaptation,
Goudarzi et al., 1997; Escudero et al., 2003; Laloë et al., by linking genetic variability to environmental data
2007), phylogeography (Hanotte et al., 2002; Matsuoka (Barker et al., 1986; Angers et al., 1999), whereas the
et al., 2002; Ciofi et al., 2006), landscape genetics (Angers impossibility of formulating an explicit model of
et al., 1999; McRae et al., 2005) and the identification of adaptation would make Bayesian clustering methods
inapplicable, in most cases. Lastly, multivariate methods
Correspondence: Dr T Jombart, UMR CNRS 5558–LBBE, ‘Biométrie et have been developed and used extensively for more than
Biologie Évolutive’, UCB Lyon 1—Bât. Grégor Mendel, 43 bd du 11 a century in various fields, such as psychometry and
novembre 1918, 69622 VILLEURBANNE cedex, France.
E-mail: jombart@biomserv.univ-lyon1.fr ecology (Pearson, 1901). Currently, multivariate analysis
Received 25 July 2008; revised 24 November 2008; accepted 8 represents a whole, rich and diversified area of statistics
December 2008; published online 21 January 2009 offering a wide choice of methods, each with its own
 Multivariate analysis of genetic markers
T Jombart et al

331
properties (Takeuchi et al., 1982; Jambu, 1991; Legendre Central to the analysis of a dataset of n objects and p
and Legendre, 1998). descriptors is the question of whether we seek a
The unfortunate consequence of this diversity of description of the relationships among the objects or
methods is that multivariate analyses are often used as among the descriptors. When analysing genetic markers,
black boxes when applied to genetic markers, leading to the main interest is in finding relationships among
frequent mistakes that sometimes question the results of objects (genotypes or populations) using p alleles. In this
an entire study. In fact, it can be difficult to know which case, data are seen as a cloud of n points embedded
method can be efficiently applied to extract information inside a p-dimensional space, where each dimension is
from genetic markers, what precautions need to be taken defined by an allele. Inside this space, inertia measures
and how the results should be interpreted. Moreover, the dispersion of n points with respect to a given
there is no doubt that multivariate analysis has been distance: this measurement of variability is used as a
under-utilized and has much more to offer to the study criterion that is optimized by the analysis. The directions
of the genetic variability. The purpose of this paper is to inside this space reflecting the highest ‘variability’ (that
critically review the use of ordination in reduced space to is, with maximum inertia) among objects are the principal
infer biological structures from genetic markers. axes, also referred to as the factors of the analysis. By
First, we attempt to clarify the rationale for these extension, a plane formed by two principal axes is often
methods and provide an overview of their current called a factorial plane. Each principal axis is defined by p
application to genetic markers. Frequent mistakes re- coordinates inside the p-dimensional space, representing
garding the utilization of these methods are then the loadings of the alleles. The principal axes are
detailed, and guidelines are provided to avoid these orthonormal (that is, perpendicular and with length
pitfalls. The following section focuses on some particula- one), and can therefore be used as a new basis to
rities of genetic markers that should be taken into represent the n objects. The set of coordinates of the
account to improve their multivariate analysis. The rest objects in this new basis are the principal components, but
of this review covers the use of multivariate analyses to the terms scores (of objects) and synthetic variables are also
tackle specific questions of interest, such as the coherence commonly used. Each principal component is associated
of the information of different genetic markers, linkage of with an eigenvalue that quantifies the amount of inertia
genetic markers to other types of data, and the study of contained in the component. Eigenvalues can also be
spatial genetic patterns. We conclude by examining some expressed as proportions of the total inertia of the
promising perspectives offered by these approaches to analysis to indicate what fraction of the entire genetic
answer challenging questions in various fields, such as variability is represented by the corresponding principal
conservation, spatial genetics and molecular ecology. components. The plot of the eigenvalues sorted in
decreasing order (the screeplot) is the basic tool used to
choose which principal components to interpret: it
Multivariate analysis of genetic markers describes how the total inertia is distributed across the
principal axes. The basic idea is that a boundary between
Rationale of multivariate analysis true structure and random noise would be indicated by a
Throughout this paper, the terms ‘ordination in reduced sharp decay between two successive eigenvalues. How-
space’ and ‘multivariate analysis’ are used interchange- ever, this is a simplistic view, and such a boundary rarely
ably. However, the first term is certainly more accurate exists in practice: the screeplot merely provides insight
than the second because ordinations in reduced space about which component likely contains interesting
represent a particular class of multivariate methods, structures, and which does not. Hence, the screeplot
another being, for instance, hierarchical clustering. The and the proportions of inertia associated with the
purpose of these methods is to summarize a strongly principal components are two complementary tools,
multivariate dataset into a small set of uncorrelated respectively indicating the genetic structures to be
synthetic variables. In other words, ordinations in reduced retained and their magnitude. The last criterion for
space aim to provide a simplified, yet meaningful, interpreting principal components is that of the biologi-
picture of complex information that is impossible to cal meaning, and is sometimes more useful than
perceive. This task implies a necessary loss of informa- statistical criteria. In some cases, the first principal
tion, and the crucial point in all these methods components (associated to large inertia) may indicate a
is to define a criterion that is optimized by the trivial structuring, and provide little biological insight.
synthetic variables seeked. For instance, in principal Conversely, principal components associated to smaller
component analysis (Pearson, 1901; Takeuchi et al., 1982, eigenvalues might contain biologically relevant informa-
pp 185–224), synthetic variables best preserve the tion; the interpretation of such components should not be
variance among observations, whereas the w2 distances discarded on the basis of a small inertia.
are preserved in the correspondence analysis (Greenacre, If multivariate analyses are unified by a single
1966). Below, we introduce general concepts required to algorithm, the core difficulty is in choosing the method
describe multivariate analyses with accuracy. that best matches the nature of the data and the questions
As formalized by the duality diagram framework asked. Because of the variety of questions and data,
(Escoufier, 1987; Dray and Dufour, 2007), most multi- numerous ordinations in reduced space are used to
variate analyses are particular cases of a general analyse genetic markers.
algorithm, and can be described using a small set of
concepts. The terminology we employ encompasses the Applications to genetic markers
most common terms, which can be found in reference Multivariate analyses are natural tools to extract
textbooks (for examples, Takeuchi et al., 1982; Jambu, biological structures from genetic markers, as these
1991; Legendre and Legendre, 1998; Lebart et al., 2004). data typically contain large numbers of genotypes or

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

332
Table 1 Multivariate analyses applied to genetic markers
Method Criterion Application Data

Principal component analysis (PCA) Variance (same as squared Euclidean Cavalli-Sforza (1966) Allozymes
distances)
Principal coordinates analysis (PCoA) Any Euclidean distance Sanchez-Mazas and Langaney (1988) Allozymes
Non-metric dimensional scaling Ordering of objects Lessa (1990) Roger’s and
(NMDS) Nei’s distances
Correspondence analysis (CA) w2 distance She et al. (1987) Allozymes
Discriminant analysis (DA) Variance between groups/total Smouse et al. (1982) Allozymes
variance
Constant-row total multiple Correlation ratio Guinand (1996) Allozymes
correspondence analysis (CRT-MCA)
Factor analysis (FA) ’Common effect’ in allele frequencies Taylor and Mitton (1974) Allozymes
Canonical correspondence analysis w2 distances in predicted data Angers et al. (1999) Microsatellites
(CCA)
Redundancy analysis (RDA) Variance of predicted data Kölliker et al. (2008) AFLP and SSR
Canonical correlation analysis (CCorA) Squared correlation between Johnson and Schaffer (1973) Allozymes
pairs of scores
Co-inertia analysis (COA) Squared covariance between Jarraud et al. (2002) AFLP
pairs of scores
Multiple co-inertia analysis (MCOA) Squared covariance between a Laloë et al. (2007) Microsatellites
set of scores
Spatial principal component analysis Product of variance and spatial Jombart et al. (2008) Microsatellites
(sPCA) autocorrelation
Abbreviations: AFLP, amplified fragment length polymorphism; SSR, single sequence repeats.
Each method is indicated by its most frequent name and abbreviation. The ‘criterion’ is the quantity optimized by the principal components
of the method. The ‘application’ column gives the reference of an early and representative publication using the method to analyse genetic
markers.

populations described by hundreds of alleles (in terms respect, NMDS can be thought of as a non-linear form of
of absolute or relative frequencies). A summary of the PCoA (Lessa, 1990). It is noteworthy that unlike other
application of these methods to genetic markers is multivariate analyses, the NMDS solution is not analy-
provided in Table 1. tical: an iterative algorithm aims at finding a good
Ordinations in reduced space are primarily used to solution, but does not guarantee that this solution is the
find a few principal components that reflect as much of best. As an alternative to PCA of allele frequencies and
the genetic variability as possible. PCA was first PCoA (or NMDS) of genetic distances, correspondence
employed to infer population structuring (Cavalli-Sforza, analysis (CA, Greenacre, 1966) can be used to analyse
1966) and spatial genetic structuring (Menozzi et al., a table of allele counts per population (She et al., 1987;
1978; Bertranpetit and Cavalli-Sforza, 1991; Cavalli- Li et al., 2002). The last multivariate analysis commonly
Sforza et al., 1993) in humans. PCA was also used early applied to genetic markers is discriminant analysis (DA,
to infer adaptations from allozyme frequencies, by Lachenbruch and Goldstein, 1979). DA is not a fully
testing the correlations between principal components exploratory approach, in that groups of genotypes must
of genetic data and principal components of a PCA of be known in advance. However, it can be used to achieve
environmental variables (Johnson et al., 1969). In disease the best discrimination between groups inside a reduced
studies, regression onto the principal components of the space, to test for genetic differentiation, and for assign-
PCA has been recently proposed to correct for popula- ment purposes (Smouse et al., 1982; Beharav and Nevo,
tion stratification (Price et al., 2006). Another method 2003).
commonly used to infer genetic structuring among Other methods have remained somewhat unnoticed,
genotypes or populations is principal coordinates analy- such as constant-row multiple correspondence analysis
sis (PCoA, Gower, 1966; Sanchez-Mazas and Langaney, (CRT-MCA, Guinand, 1996; Guinand et al., 1996), factor
1988; Warnes, 2003). Although PCA preserves the analysis (FA, Taylor and Mitton, 1974; Mulley et al., 1979)
canonical Euclidean distance among the studied entities, and distance-based redundancy analysis (db-RDA,
PCoA can be employed to summarize any Euclidean Legendre and Anderson, 1999; Geffen et al., 2004).
genetic distance between genotypes or populations, but The reason for this may be historical, or could arise
does not provide a representation of the alleles. This from problems associated with using these approaches.
offers the advantage of using measures of genetic For instance, CRT-MCA aims at finding synthetic
variability that are directly related to a population variables with maximum FST, but only proposes
genetics model; for instance, PCoA has been used to an approximate solution. Denoting f as a set of
summarize matrices of pairwise FST (Zhivotovsky et al., frequencies of an allele for q populations, f̄ as the mean
2003) and of Roger’s distance (Baker and Moeed, 1987). frequency computed across populations and var( f ) as
Non-metric dimensional scaling (NMDS, Cox and Cox, the variance between populations of f, FST is defined
2001) has also been employed to analyse matrices of as var( f )/f̄ (1
f̄ ), where f̄ (1
f̄ ) is the theoretical
genetic distances (Baker and Moeed, 1987; Lessa, 1990). variance of f (Weir, 1996, p 166). Unfortunately,
However, NMDS differs from PCoA in that it attempts to the quantity optimized by CRT-MCA is var( f )/s2f ,
preserve the ordering of objects based on their genetic where P s2f is the empirical variance of f
q
distance rather than their genetic distance per se; in this (s2f ¼ 1q i¼1 ðfi
fÞ2 Þ:While for arbitrarily large samples

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

333
s2f converges towards f̄ (1
f̄ ), these quantities differ in populations of common minas (Acridotheres tristis). Lack
practice, and the principal components yielded by CRT- of accuracy in the description of the method always
MCA do not optimize FST. A possible cause for the minimal complicates interpretation of the results, and sometimes
use of FA is that it was introduced to correlate patterns in brings their validity into question. For instance, some
allele frequencies with environmental variables (Taylor and papers show principal components of a PCA that were
Mitton, 1974), which is not the purpose of this method. In clearly not centred (their range of variation did not
fact, FA estimates a model in which allele frequencies are include zero), which indicates an error in the computa-
expressed as a sum of two components: a part common to tions of the analysis and invalidates the results (for
every allele and a residual part representing allele-specific instance, MacHugh et al., 1997, 1998; Pariset et al., 2003).
effects (Seal, 1966, pp 153–180). Lastly, it is not clear why Moreover, it is difficult to ascertain precisely where the
db-RDA has not been applied more often to genetic problem came from, as the software used for the
markers, but this could simply be due to its recent computation was not mentioned in these publications.
application (Geffen et al., 2004).
Although multivariate analyses can be efficiently used
Making graphics
to extract information from genetic markers, choosing a
Another classical problem lies in the graphical display of
method appropriate to the data and the question being
results. As mentioned previously, the screeplot is the
asked is sometimes difficult. As a matter of fact, a
basic tool used to assess which principal components
number of mistakes occur quite frequently in such
should be interpreted, but it is most often omitted in
applications. In the following, we point out the major
publications. The amount of inertia associated with each
pitfalls, as well as strategies to avoid them.
principal component is often indicated, but this informa-
tion is complementary to the screeplot and cannot be
Misuses, misinterpretations and specific used as a substitute. For instance, in their study of the
genetic differentiation among different yak (Poephagus
issues grunniens) populations, Xuebin et al. (2005) presented
a scatterplot of PCA displaying 80% of the whole
Ensuring reproducibility
variability, but this scatterplot was merely uninformative
A first concern in data analysis is to ensure reproduci-
in terms of genetic differentiation. Conversely, two
bility, or at least to provide all the elements required to
principal components of PCA containing less than 10%
evaluate the relevance of the results. Unfortunately, the
of total inertia provided insights about the phylogeny of
literature regularly provides examples of studies in
different maize subspecies in Matsuoka et al. (2002).
which it is almost impossible to know which analyses
When used alone, the amount of inertia can therefore be
were actually performed.
a misleading criterion for choosing the principal compo-
The first problem lies in the absence of an accurate
nents to interpret (see ‘Interpreting genetic structures’).
description of the method used: reference articles are
Another widespread custom is the use of 3-dimen-
rarely cited, and abbreviations sometimes do not match
sional scatterplots (van Pijlen et al., 1995; Xuebin et al.,
the name of the method. For instance, ‘PCA’ is used to
2005). Although these representations add a fancy touch
refer to principal coordinates analysis (PCoA) in Pariset
to multivariate analyses, they also have the unfortunate
et al. (2003). Such confusion adds to the ambiguities that
effect of sacrificing the mathematical properties of an
already exist between some methods, such as those
analysis, and thus its interpretability. By definition,
between PCoA and NMDS. PCoA is also sometimes
principal axes and the associated principal components
called ‘metric dimensional scaling’ (MDS), whereas
provide the best possible planar representation of the
NMDS is indifferently abbreviated MDS or NMDS
data. If three principal components are retained, their
(Legendre and Legendre, 1998). This is all the more
representation requires two factorial planes, with one
confusing since PCoA is routinely used to initialize the
axis being redundant. Scatterplots in three dimensions
algorithm of NMDS (Baker and Moeed, 1987). Papers
are ultimately always viewed on a screen or on a sheet of
demonstrating an ambiguity between PCoA and NMDS
paper, and are thus re-projections of three principal
are not uncommon (for example, Preziosi and Fairbairn,
components in two dimensions. The obtained represen-
1992; Zhivotovsky et al., 2003).
tations are necessarily worse than the true representation
Although required, providing a correct reference to a
of principal components because they no longer have
method is usually not sufficient. Some methods exist in
maximum inertia nor orthogonality. Hence, 3-dimen-
different variants, according to the initial transform-
sional visualization should be restricted to interactive
ations of the data. This is particularly true for PCA:
data analysis (where it can be useful), and is better
although centring (subtracting the mean allele frequency
avoided in publications.
from all observations) is always achieved, scaling of the
Apart from these pitfalls common to every multi-
alleles (dividing each observation by allelewise values) is
variate analyses, some specific issues also arise when
optional and can be performed in several ways. Scaling
certain methods are applied to genetic markers.
can drastically change the results of a PCA, but is rarely
disclosed (for example, Mitton, 1978; MacHugh et al.,
1998; Grivet et al., 2008). In PCoA and NMDS, the genetic Some specific issues
distance employed should always be specified, and in A first particular issue concerns the use of CA. This
the case of NMDS, how the algorithm was initialized method is appropriate for the analysis of a contingency
should be indicated. An example of such an application table, that is, a matrix of positive integers (Greenacre,
can be found in Baker and Moeed (1987), who used an 1966), and is thus appropriate for the analysis of a table
NMDS initialized by a PCoA of Roger’s distances of of allele counts. A good example of such an application is
allozyme data to explore the genetic variation among provided by She et al. (1987), who used a CA of allozyme

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

334
data to investigate the genetic differentiation between subset of best discriminating alleles (Fahima et al., 1999;
populations of teleost fishes. Interestingly, this study also Beharav and Nevo, 2003). However, investigations
showed that ‘correspondences’ highlighted by CA can should be carried out to assess whether a particular
reflect linkage disequilibrium existing between alleles. In variable selection procedure is preferable to the others in
some cases, CA has been used for allele (relative) the case of allele frequencies.
frequencies (Li et al., 2002), which has been proven to
significantly alter the results of the analysis (Perrière and
Thioulouse, 2002). In such a case, it seems much more Interpreting genetic structures
appropriate to use PCA or PCoA. However, even when A major concern in multivariate analysis of genetic
CA is correctly used, one should be aware that scarce markers lies in interpreting the results. This issue can be
descriptors are given a stronger weight by the w2 illustrated by examining one case of misinterpretation,
distance, which is optimized by the analysis (Legendre raising the question of which result of a multivariate
and Legendre, 1998, p 285). The typical consequence is analysis could be interpreted as genetic structuring.
that a population possessing a rare allele would appear In the controversy regarding the relevance of definiting
as an outlier in CA components. Simple simulations human races based on genetic information (Lewontin, 1978;
show that such an artifactual pattern arises even when Mitton, 1978; Powell and Taylor, 1978), Mitton (1978)
studying groups of genotypes are randomly chosen from argued that genetic differentiation between ‘human
the same population (results not shown). A way to avoid races’ was important because they clearly appeared as
this problem is to remove rare alleles from the data prior distinct groups on the factorial map of a PCA. This
to the CA, although this solution requires some misinterpretation of the results is related to the common
investigations regarding which frequency should be mistake of not displaying the screeplot of the analysis
considered as ‘rare’ from the point of view of the CA. along with the values of inertia associated with principal
A second specific issue occurs in DA. This method components. Ordinations in reduced space do not
finds principal components maximizing the variance summarize the essential part of the genetic variability:
between populations while keeping the variance inside they attempt to show as much genetic variability as
populations constant, assuring optimal discrimination possible in a few axes, which is different. Mitton (1978)
between the populations (Krzanowski and Marriott, showed that ‘racial’ groups were well separated on the
1995, pp 1–56). However, this method involves computa- factorial plane and that two principal components were
tion of the Mahalanobis metric (Beharav and Nevo, sufficient to assign each population to a given group.
2003), which is the inverse of the matrix of covariances However, this did not contradict the well-acknowledged
between alleles. For this inverse to exist, the covariance fact that the genetic variability within ‘races’ is much
matrix must be of full rank, that is, of rank p if there are p larger than between ‘races’ (Edwards, 2003), as sug-
alleles (Harville, 1997, p 80). This is never the case for gested by the author. For example, it would be possible
allele frequencies: each marker spans a space of at most to perfectly discriminate two populations using only one
one dimension less than the number of its alleles because allele, but this allele may represent only 1% of the
any frequency is entirely defined by all the others. That variability of a dataset containing 100 alleles. This point
is, if there are k markers, the rank of the covariance was discussed by (Edwards, 2003), who emphasized the
matrix is at most min(p
k, n). Thus, the discriminant fundamental difference between being able to assign
analysis can only be performed on a matrix of allele genotypes to taxonomic groups, and observing larger
frequencies after removing a given number of alleles, genetic variability between these taxonomic groups.
and assuring that there are more objects (genotypes or We can ask what criterion the principal components of
populations) than alleles. In fact, the number of objects n an analysis should meet to be considered as true genetic
should be consequently larger than the number of alleles structuring. The relative amount of inertia cannot be
p: (Williams and Titus, 1988) reported that n should be at used as a single criterion, because it depends directly on
least three times larger than p for DA to yield reliable the number of alleles considered. As stressed previously,
results. Multicollinearity can also exist among alleles the screeplot can be used to assess which principal
(that is, when alleles are correlated), especially when components likely contain interesting structures. Re-
linkage disequilibrium occurs. In these cases, the cently, Patterson et al. (2006) tested the significance of the
Mahalanobis metric is said to be ill-conditioned, result- eigenvalues from a PCA of genetic markers to infer
ing in numerical instability. As a result, principal axes population stratification. Another testing approach to
and principal components of DA cannot be computed select interpretable principal components of a PCA has
with accuracy, and small changes in allele frequencies been proposed by (Dray, 2008), and could also be used to
induce large changes in the results (Seber, 1977, pp 319– identify significant genetic structures. It is noteworthy
322). As a consequence, the alleles used in DA should be that both approaches are reserved to PCA (Patterson
carefully selected before performing the analysis. An et al., 2006; Dray, 2008), and it would be valuable to
empirical approach consists of retaining only the most extend these tests to other multivariate methods.
frequent allele of each locus (Sagnard et al., 2002), but this Another way of assessing relevant genetic structures
does not ensure that the subset of alleles obtained is emerging from an ordination method is to quantify the
optimal with regard to discrimination. In fact, it does not amount of genetic differentiation contained in the
even ensure that the multicollinearity problem is solved, principal components. The main difficulty is then
as linkage disequilibrium can still exist between the most identifying clusters of genotypes from the principal
frequent alleles. One can preferentially use statistical components retained. This can be achieved using a given
approaches that are especially devoted to the selection clustering algorithm (Legendre and Legendre, 1998,
of variables in DA (Lachenbruch and Goldstein, 1979), pp 303–381), such as the unweighted arithmetic average
where such approaches proved useful for selecting a clustering (UPGMA, Rohlf, 1963). It is then possible to

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

335
measure the amount of genetic differentiation between (Dolédec and Chessel, 1987), which maximizes the
the obtained clusters of genotypes using classical variance between populations, would yield principal
approaches like FST. Note that the obtained statistics components with maximum FST if performed on allele
can only be used to quantify genetic differentiation, but frequencies
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffifficentred
ffi to a mean of zero and scaled
not to test it, because the principal components are by by fj ð1
fj Þ . Even though between-class PCA has
definition, optimized with regard to some measurement only recently been applied to genetic markers by Parisod
of genetic differentiation. To conclude this point, the and Christin (2008, presented as ‘inter-class PCA’) and
identification of interpretable structures is a major Jombart (2008), this method seems promising for
question in multivariate analysis, and is of particular investigating genetic differentiation between groups of
interest when seeking genetic structures from molecular genotypes.
markers.
As we have seen, the application of multivariate Compositional data
analysis to genetic markers can be improved by avoiding The principal particularity of allele frequencies may be
a number of pitfalls. However, further improvements can that they are sets of compositional data, that is, data with
be gained by adapting multivariate methods to several a constant sum for each population and locus. This
particularities of genetic markers. feature induces non-independence between allele fre-
quencies inside each locus (a frequency can always be
deduced from all the others), and has several conse-
Respecting the very nature of data quences on ordinations in reduced space. Developments
in the multivariate analysis of compositional data were
Scaling in PCA led by the work of Aitchison (Aitchison, 1983, 1999, 2003;
In many cases, genetic markers are analysed as allele Aitchison and Greenacre, 2002), but remained mostly
frequencies, which are subjected to a PCoA or a PCA. ignored in genetics, apart from a few exceptions
PCoA is usually well suited to genetic markers because (Romano et al., 2003; Reyment, 2005). As stressed before,
several genetic distances can be used to summarize the allele frequencies at a given locus are not independent, as
genetic variability. In this case, it is necessary to use a one can be entirely deduced from the others. Populations
Euclidean distance like Roger’s (Weir, 1996, p 197), so described by pj alleles at the jth locus are not embedded
that genetic relationships among entities can be fully inside a pj-dimensional space, but are instead inside a
represented in a plane and to choose a distance whose space whose maximum dimensionality is (pj
1), known
underlying model best matches the data (see for instance as a simplex space (Aitchison, 2003, pp 24–28). A variety
Weir, 1996, pp 190–198). In the case of PCA, attention of problems can occur when directly computing an
must be devoted to the transformations of data: if ordination in reduced space in the simplex space (or in a
centring of allele frequencies is almost mandatory, the set of simplex spaces in the case of several loci), like the
scaling of allele frequencies can be discussed. The impossibility of identifying structures that are intrinsi-
general reason for scaling is to compensate for trivial cally non-linear and the numerical instability of principal
differences that occur in the variance of the descriptors, components. The solutions proposed to account for these
for instance, when descriptors are expressed in different problems rely on transforming frequencies (mostly using
units. A reason for not scaling allele frequencies is that logarithms) and performing a classical analysis like PCA
doing so is not necessary (scales of variation are of the obtained data. Reyment (2005) showed that the
inherently the same for every allele), and could mask results of PCA could be strikingly improved by such
differences in the genetic variability contained by practices, even when considering a simple log trans-
informative and non-informative markers, ultimately formation of the data. Henceforth, these approaches
hiding structures in the data. Nonetheless, one good should be considered when analysing allele frequencies.
argument for scaling allele frequencies would be to
compensate for differences in variance among alleles due
to their underlying binomial nature: the theoretical
Diversity inside the diversity
variance associated with the jth allele frequency, fj A portion of the literature in conservation biology
(j ¼ 1,yp where p is the total number of alleles), is stresses the idea that different genetic markers can
proportional to fj(1
fj). The result is that the variance of provide different information about the genetic diversity
an allele frequency is expected to be ‘naturally’ larger for of a set of populations (Moazami-Goudarzi and Laloë,
frequencies close to 0.5, and smaller for frequencies close 2002). In fact, genetic markers are usually taken as a
to 0 or 1. The PCA seeking linear combinations of alleles whole to seek a global, common typology of individuals
with maximum variance, alleles with frequencies closer or populations, without trying to assess if such a
to 0.5 would be favoured by the analysis, although common typology exists. There are, however, good
not necessarily reflecting a genetic structure. One
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ffi reasons for this typology not to occur, the first being
way to correct this is to divide fj by fj ð1
fj Þ , that selection can affect different loci in different ways. If
as has been previously proposed (Cavalli-Sforza et al., this is obvious for selected markers like allozymes, it can
1994, pp 41–42). Mulley et al. (1979) used a related also be true for supposedly neutral markers that are
standardization of allele frequencies, which does not physically linked to selected regions of the genome.
amount to unit theoretical variance, but accounts for Interestingly, the first studies linking the genetic varia-
the number of genotypes used to compute frequencies bility in allozymes to environmental features analysed
in each population. Interestingly enough, the variance each locus separately by PCA (Johnson et al., 1969;
between populations of
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ffi the allele frequency Johnson and Schaffer, 1973).
standardized by fj ð1
fj Þ is exactly the classical FST To tackle the question of the typological coherence of
(Weir, 1996, p 166). Therefore, the between-class PCA genetic markers, the locus must be considered as the unit

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

336
of analysis. In this perspective, if there are K markers, methods presented below (their purpose is to investigate
K analyses should be performed and compared. A class the genetic differentiation between groups of genotypes),
of multivariate analyses, called the K-table methods DA and between-class PCA are not presented in this
(Dray et al., 2007), is devoted to this particular task. Such section.
methods were introduced in genetics by Laloë et al.
(2007), who used multiple co-inertia analysis (Chessel Asymmetric methods
and Hanafi, 1996) to compare the typological information The first type of method is formed by constrained
provided by different microsatellites. This study showed ordinations, which are devoted to investigating the
that microsatellites could provide different pictures of variability in one dataset that can be explained by
the genetic diversity among populations: whereas some another dataset. This is achieved by a multivariate
microsatellites reveal the entire genetic structure, some regression of a ‘response’ dataset onto an ‘explanatory’
perceive only particular aspects of the genetic diversity dataset (Ter Braak, 1986). These methods are thus
and others are simply not informative in terms of genetic asymmetric, in that the variability in one dataset is
differentiation patterns. The typological value of a explained by another. There are two main techniques in
marker can be used to quantify the extent to which this this context: redundancy analysis (RDA, Rao, 1964),
marker contributes to displaying a particular genetic which is a constrained version of PCA, and canonical
structure (Laloë et al., 2007). The application of K-table correspondence analysis (CCA, Ter Braak, 1986), which
approaches to genetic markers was further developed by is based on CA. RDA and CCA therefore inherit their
(Pavoine and Bailly, 2007), who introduced other properties from PCA and CA: RDA can be used for allele
K-table methods coupled with a multivariate analysis frequencies, whereas CCA is more appropriate to
of biodiversity (Pavoine et al., 2004). Their results analysis of tables of allele counts. Both RDA (Kölliker
confirmed the fact that summing the information coming et al., 2008) and CCA (Angers et al., 1999) have proven
from different genetic markers, as is usually performed useful in population genetics, mostly to investigate the
for ordinations in reduced space, does not always portion of the genetic variability that can be explained by
provide the most accurate picture of biodiversity. Note a set of environmental variables. For instance, in Angers
that if K-table methods can suggest that loci experience et al. (1999), the CCA revealed that the genetic diversity
different selective pressures, they cannot be used as a among a set of brook charr populations (Salvelinus
direct test for these differences. In fact, K-table ap- frontalis) was mainly driven by the structure of the
proaches are first and foremost designed to identify hydrographic network and by a few environmental
common typologies, and not discrepancies, among a set variables. Another interest of this study is that analyses
of markers. were applied to two different levels, to study the effects
If K-table methods are more complex tools than of hydrographic and environmental features on the
single-table analyses, their use in genetics should be genetic diversity inside, and between populations.
considered with attention. Note that the linkage of Like discriminant analysis, RDA and CCA involve
multilocus genetic information to environmental computation of the Mahalanobis metric which is, in this
features like in Johnson et al. (1969) still raises challen- case, the matrix of covariances between explanatory
ging questions in terms of data analysis: How can we variables (Legendre and Legendre, 1998). These analyses
describe the genetic-environment relationships at therefore require that the number of explanatory variables
several loci? What are the different patterns of adapta- (for instance, environmental variables) be fairly lower than
tion among loci? the number of studied objects (genotypes or populations)
to be computable. Following the previously cited study of
Williams and Titus (1988) concerning DA, we can
Linking genetic markers to other data recommend that the number of objects should be at least
One of the greatest applications of ordinations in three times larger than the number of explanatory
reduced space is in the linkage of genetic markers to variables. RDA and CCA also demand that the explanatory
other types of data (Johnson et al., 1969; Taylor and variables are reasonably uncorrelated to achieve numerical
Mitton, 1974; Mulley et al., 1979; Barker et al., 1986; stability and interpretability of the results. As a rule of
Jarraud et al., 2002). This is typically the case in the study thumb, we suggest avoidance of correlations greater than
of genotype-environment relationships, where multi- 0.7, so that no more than one half of the variability of any
variate methods can be used to investigate correlations predictor could p beffiffiffiffiffiffi
explained
ffi by another predictor (that is,
between genetic data and environmental features (John- R2 o0:5 , jrjo 0:5 ’ 0:7 ). Note that genetic markers
son et al., 1969; Mulley et al., 1979). Another application could also be used as explanatory variables, for example,
of such an approach is to relate genetic information to with an ‘explained’ dataset of phenotypic traits. In such
phenotypic data (Jarraud et al., 2002). Note that when cases, the dimension of the genetic information should be
patterns of selection are being investigated, the genetic reduced, either by applying a standard variable selection
diversity should be inferred from non-neutral rather than procedure (for example, forward selection) to the allele
neutral markers. Various methods are available for frequencies, or by reducing the genetic data to a few
coupling two different kinds of information, some of principal components using PCA or PCoA.
which have been introduced into population genetics. When the above conditions are respected, constrained
These can be divided into two categories, depending on ordinations can be efficiently used to explain one kind of
whether they treat both types of information symme- variability by another. However, when the purpose of a
trically or not. Approaches like DA and between-class study is to investigate common patterns of variability in
PCA are also methods for coupling genetic markers with two datasets, or when RDA and CCA cannot be used for
a different information (some partitions of individuals). technical reasons, an alternative can be found in certain
However, because their aim is very different from the symmetric approaches.

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

337
Symmetric methods analysis of genetic markers coupled to other kinds of
Symmetric methods allow one to study the structures information (Cavalli-Sforza et al., 1994, p 41), although
common to two datasets by treating the two types of we were unable to find any applications of this technique
information similarly. They differ from constrained to genetic markers. COA finds two sets of principal axes
ordinations in the same way that linear regression differs (one for each dataset), such that the pairs of principal
from correlation. Symmetric approaches include canoni- components have a maximum squared covariance
cal correlation analysis (CCorA, Hotelling, 1936; Takeu- (that is, co-inertia). This criterion is particularly interest-
chi et al., 1982, pp 225–280) and co-inertia analysis (COA, ing as it amounts to maximizing the product of the
Dolédec and Chessel, 1994; Dray et al., 2003a). CCorA variances of each principal component and their squared
was introduced by Johnson and Schaffer (1973) to correlations (because cov2(a,b) ¼ var(a)var(b)cor2(a,b)).
describe and test the correlations between allele frequen- Interestingly, the COA does not require inversion of a
cies in allozymes and a set of environmental features. covariance matrix; consequently, it does not require the
The principle of this technique is to find two sets of number of descriptors to be lower than the number of
orthogonal axes (one for each dataset), such that the objects and it is not hampered by correlations among the
obtained pairs of principal components have a maximum descriptors. Moreover, COA relies on a modification of
squared correlation (Takeuchi et al., 1982, pp 225–229). It two separate analyses, each of which can be, for instance,
is worth noting that Johnson and Schaffer (1973) were a PCA, a PCoA, or a CA. For example, Jarraud et al.
following another pioneering work (Johnson et al., 1969) (2002) employed the co-inertia between a PCoA of a
in which the same authors used correlations between genetic distance matrix derived from AFLP markers and
principal components of two PCAs (one of allele a PCA of distributions of toxin genes in several strains of
frequencies, one of environmental variables) to test Staphylococcus aureus to assess the evolution of virulence
genetic-environment relationships. A series of subse- factors with respect to the genetic background of the
quent papers provided remarkable illustrations of the strains. The COA appears to be a good alternative to
insights that CCorA can bring to the study of adaptation RDA, CCA and CCorA when these methods cannot be
(Schaffer and Johnson, 1974; McKechnie et al., 1975; applied for the reasons described above. In other cases,
Mulley et al., 1979). A nice example is provided by the COA may still be favoured whenever the squared
Mulley et al. (1979), which used the CCorA to investigate covariance criterion is more satisfying than criteria used
patterns of adaptation in populations of Drosophila by other analyses, that is, when one is interested in
buzzatii. The authors have shown that the allelic variation identifying common patterns of variation between two
observed at some allozyme loci was significantly different sources of information.
correlated to climate descriptors, which strongly sug-
gested the existence of local adaptations in these
populations. A recurrent problem in these studies is that Spatial multivariate analysis
gene-flow can act as a confounding effect when assessing Many population genetics studies in which multivariate
genetic-environment correlations. Schaffer and Johnson analyses were used involve georeferenced data. When
(1974) addressed this issue by regressing allele frequen- processes related to gene-flow are being investigated—
cies onto spatial coordinates prior to the analysis, and which may be the most common case—spatial genetic
hence removing linear spatial trends from the data. Note patterns are researched in neutral markers (Menozzi
that more efficient methods of removing spatial patterns et al., 1978; Cavalli-Sforza et al., 1993). In contrast, when
have since been developed, some of which are described non-neutral markers are used to infer patterns of
in the next section. adaptation, spatial structures induced by gene flow can
A typical problem in CCorA is that, like RDA and act as a confounding effect that would have to be
CCA, it requires to compute the Mahalanobis metric of removed (Schaffer and Johnson, 1974). As noted by
both datasets: it cannot be used when there are more Mulley et al. (1979), the drawback of this strategy is that
descriptors than studied objects and it requires descrip- ‘if environmental factors with selective effects are strongly
tors to be uncorrelated to yield interpretable results. In correlated to geographic location, adjustment for location may
some of these cases, a CCorA can still be performed after remove a major fraction of the selective effects’. In such a case,
selecting a small subset of uncorrelated variables (for it would be worthwhile to compare the selective effects
example, Mulley et al., 1979). A common criticism of detected with and without removing the effects of spatial
CCorA is that pairs of principal components with patterns. Spatial information can be used in multivariate
maximum squared correlation could have a very small analysis of genetic markers, to investigate the part of the
variance, and therefore have in general no real biological genetic variability that is or is not spatially structured.
meaning (Taylor and Mitton, 1974). Taylor and Mitton Unfortunately, the methods commonly used to inves-
(1974) suggested that a symmetric analysis should yield tigate spatial genetic patterns almost never take
pairs of principal components reflecting both a fair spatial information into account explicitly, that is, they
amount of variance and be correlated with each other, do not incorporate spatial information as a component of
that is, reflecting common parts of the variability in the the criterion optimized by the analysis (Jombart et al.,
two datasets. This is the definition of a method 2008). This contrasts with other methodological
developed later in ecology; the co-inertia analysis frameworks such as analysis of molecular variance
(Dolédec and Chessel, 1994; Dray et al., 2003a). (Excoffier et al., 1992) or Bayesian clustering (Pritchard
COA has been imported into genetics to relate the et al., 2000), in which spatially explicit methods are used
genetic variability of several bacterial strains to the (respectively, Dupanloup et al., 2002; François et al., 2006).
expression of toxin genes (Jarraud et al., 2002). It is worth However, spatial ordinations exist and are widely used
noting that COA is closely related to Procustean analysis in other domains, the closest to genetics being ecology. It
(Dray et al., 2003b), which has been proposed for the is, therefore, not surprising that spatial ordinations were

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

338
first proposed to analyse genetic markers in vegetation an appealing choice: in addition to allowing exact
sciences (Escudero et al., 2003) and landscape genetics reproducibility, it provides an interface between a large
(Grivet et al., 2008). number of implemented multivariate methods (Chessel
Recently, Grivet et al. (2008) used the canonical et al., 2004; Dray et al., 2007) and genetic marker data
trend surface analysis (Wartenberg, 1985) to detect (Jombart, 2008), in addition to supporting the usual
spatial patterns using microsatellite markers. This population genetics tools (Warnes, 2003; Goudet, 2005).
approach relies on performing a CCorA to identify From a more theoretical point of view, it seems important
correlations between genetic and spatial data. Grivet to further investigate the relationships between multi-
et al. (2008) used polynomials of spatial coordinates as variate methods and genetic models. A step in this
spatial predictors, although this approach was criticized direction has been made by Patterson et al. (2006), who
in ecology (Borcard and Legendre, 2002; Dray et al., applied recent developments in statistics (Soshnikov and
2006), mainly because the obtained variables are gen- Fyodorov, 2005) to infer the number of populations in a
erally correlated and can only model broad-scale set of genotypes and define a threshold for genetic
patterns. Other spatial predictors, Moran’s eigenvectors, structuring to be detectable by PCA.
are now used in ecology (Dray et al., 2006; Griffith and More generally, several multivariate analyses devel-
Peres-Neto, 2006). Contrary to polynomials of spatial oped in other disciplines can be adapted to search
coordinates, these spatial predictors are uncorrelated, biological structures within genetic markers. This is
and can model spatial patterns on a wide range of scales. clearly the case in spatial genetics, where constrained
To reveal spatial genetic patterns, Moran’s eigenvectors ordinations based on Moran’s eigenvectors (Dray et al.,
can be used as explanatory variables in a CCA or an 2006) could be used to investigate or correct for
RDA of genetic markers. In studies in which spatial spatial genetic structures. It is also true for K-table
structures need to be removed to infer adaptations, methods, which were only recently introduced into
Moran’s eigenvectors could also be used as covariables population genetics (Laloë et al., 2007; Pavoine and
in partial RDA or partial CCA (Legendre and Legendre, Bailly, 2007), and open appealing perspectives for the
1998, pp 769–779). study of the genetic diversity. These methods can also be
To our knowledge, the only spatial ordination devel- used to investigate common patterns of variation
oped within the genetic framework is the spatial inferred from genetic markers and other sources of
principal component analysis (sPCA, Jombart et al., information, like biological traits and environmental
2008). This method relies on a modification of PCA such features. As noted by Patterson et al. (2006), multivariate
that not only the variance of the principal components, analysis can analyse larger datasets than other usual
but also their spatial autocorrelation, is optimized. approaches such as Bayesian clustering, and thus
Jombart et al. (2008) identified various kinds of represents a relevant approach to extracting information
spatial structuring that can arise in genetic data, and from huge datasets produced by the detailed mapping of
showed that sPCA can be efficiently used to reveal genetic variation for a large number of genotypes. This is
these patterns. In particular, a comparison between the case, for instance, with the ‘1000 Genomes’ project
PCA and sPCA demonstrated that sPCA should be (http://www.1000genomes.org/), which aims at se-
preferred to PCA whenever spatial genetic patterns quencing one thousand human genotypes to provide
are researched. Note that a similar approach was high-resolution information that is directly valuable for
developed in the vegetation sciences by Dray et al. disease studies. Promisingly, a wide range of questions
(2008), who proposed a spatial version of CA. Although are raised by or through genetic markers, some of which
sPCA is devoted to investigating spatial genetic patterns can currently be solved by existing methods. Some of
in allele frequencies, the approach of Dray et al. (2008) these questions will undoubtedly require specific devel-
could be used to study spatial genetic patterns in allele opments in which multivariate models will have to
counts. closely match the genetic concerns, which makes the
multivariate analysis of genetic markers a whole area of
research in biometry.
Perspectives and conclusion
We reviewed how a multivariate analysis can be used to
extract biological information from genetic markers. The Acknowledgements
large diversity of existing multivariate methods allow to We are very grateful to Christian Biémont, F Stephen
investigate a wide variety of genetic structures, which Dobson, Gilles Yoccoz and Sébastien Devillard for
depend on the nature of data as well as on the question providing constructive comments on an earlier version
being asked. One important observation emerging from of the article. We also address many thanks to F Stephen
this review is that application of multivariate methods to Dobson for improving the English of this paper.
genetic markers could sometimes benefit from more
rigorous practices. Methods should always be referred to
clearly and with a distinction between the method itself References
and its implementation. An accurate description of an
ordination in reduced space would include all data Aitchison J (1983). Principal component analysis of composi-
tional data. Biometrika 70: 57–65.
transformations, such as centring and scaling in PCA, the
Aitchison J (1999). Logratios and natural laws in compositional
chosen distance in PCoA and NMDS, the selection of data analysis. Math Geol 31: 563–589.
alleles in DA or algorithm initialization in NMDS. To Aitchison J, Greenacre M (2002). Biplot of compositional data.
facilitate reproducibility, free and script-based software J R Stat Soc Ser C, Appl stat 51: 375–392.
should be favoured over other software. In this context, Aitchison J (2003). The Statistical Analysis of Compositional Data.
R software (R Development Core Team, 2008) is clearly The Blackburn Press: Cladwell, New Jersey.

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

339
Angers B, Plante M, Bernatchez L (1999). Canonical correspon- Escudero A, Iriondo JM, Torres ME (2003). Spatial analysis of
dence analysis for estimating spatial and environmental genetic diversity as a tool for plant conservation. Biol Conserv
effects on microsatellite gene diversity in brook charr 113: 351–365.
(Salvelinus fontinalis). Mol Ecol 8: 1043–1053. Excoffier L, Smouse PE, Quattro JM (1992). Analysis of
Baker AJ, Moeed A (1987). Rapid genetic differentiation and molecular variance inferred from metric distances among
founder effect in colonizing populations of common mynas DNA haplotypes: application to human mitochondrial DNA
(Acridotheres tristis). Evolution 41: 525–538. restriction data. Genetics 131: 479–491.
Barker JSF, East PD, Weir BS (1986). Temporal and microgeo- Fahima T, Sun GL, Beharav A, Krugman T, Beiles A, Nevo E
graphic variation in allozyme frequencies in a natural (1999). RAPD polymorphism of wild emmer wheat
population of Drosophila buzzatii. Genetics 112: 577–611. populations, Triticum dicoccoides, in Israel. Theor Appl Genet
Beharav A, Nevo E (2003). Predictive validity of discriminant 98: 434–447.
analysis for genetic data. Genetica 119: 259–267. Falush D, Stephens M, Pritchard JK (2003). Inference of
Bertranpetit J, Cavalli-Sforza LL (1991). A genetic reconstruc- population structure using multilocus genotype data:
tion of the history of the population of the Iberian Peninsula. linked loci and correlated allele frequencies. Genetics 164:
Ann Hum Genet 55: 51–67. 1567–1587.
Borcard D, Legendre P (2002). All-scale spatial analysis of Fisher RA (1952). Statistical methods in genetics. Heredity 6:
ecological data by means of principal coordinates of 1–12.
neighbour matrices. Ecol Modell 153: 51–68. François O, Ancelet S, Guillot G (2006). Bayesian clustering
Cavalli-Sforza LL (1966). Population structure and human using hidden markov random fields in spatial population
evolution. Proc R Soc Lond Ser B 164: 362–379. genetics. Genetics 174: 805–816.
Cavalli-Sforza LL, Menozzi P, Piazza A (1993). Demic expan- Geffen E, Anderson MJ, Wayne RK (2004). Climate and habitat
sions and human evolution. Science 259: 639–646. barriers to dispersal in the highly mobile grey wolf. Mol Ecol
Cavalli-Sforza LL, Menozzi P, Piazza A (1994). The History 13: 2481–2490.
and Geography of Human Genes. Princeton University Press: Goudet J (2005). HIERFSTAT, a package for R to
Princeton. compute and test hierarchical F-statistics. Mol Ecol Notes 5:
Chessel D, Hanafi M (1996). Analyses de la co-inertie de K 184–186.
nuages de points. Revue de statistique appliquée XLIV 2: Gower JC (1966). Some distance properties of latent root and
35–60. vector methods used in multivariate analysis. Biometrika 53:
Chessel D, Dufour AB, Thioulouse J (2004). The ade4 package-I- 325–338.
one-table methods. R News 4: 5–10. Greenacre M (1966). Theory and Applications of Correspondence
Ciofi C, Wilson GA, Beheregaray LB, Marquez C, Gibbs JP, Analysis. Academic Press: London.
Tapia W et al. (2006). Phylogeographic history and gene flow Griffith DA, Peres-Neto P (2006). Spatial modeling in ecology:
among giant galápagos tortoises on southern Isabela Island. the flexibility of eigenfunction spatial analyses. Ecology 87:
Genetics 172: 1727–1744. 2603–2613.
Cox RF, Cox MAA (2001). Multidimensional Scaling. Chapman & Grivet D, Sork VL, Westfall RD, Davis FW (2008). Conserving
Hall/CRC: Bora Raton, Florida. the evolutionary potential of California valley oak (Quercus
Dolédec S, Chessel D (1987). Rythmes saisonniers et compo- lobata Née): a multivariate genetic approach to conservation
santes stationnelles en milieu aquatique. I. description d’un planning. Mol Ecol 17: 139–156.
plan d’observation complet par projection de variables. Acta Guinand B (1996). Use of a multivariate model using allele
Oecologica, Oecologia Generalis 8: 403–426. frequency distributions to analyse patterns of genetic
Dolédec S, Chessel D (1994). Co-inertia analysis: an alternative differentiation among populations. Biol J Linnean Soc 58:
method for studying species-environment relationships. 173–195.
Freshw Biol 31: 277–294. Guinand B, Bouvet Y, Brohon B (1996). Spatial aspects of genetic
Dray S, Legendre P, Peres-Neto P (2006). Spatial modelling: a differentiation of the European chub in the Rhone River
comprehensive framework for principal coordinate basin. J Fish Biol 49: 714–726.
analysis of neighbour matrices (PCNM). Ecol Modell 196: Hanotte O, Bradley DG, Ochieng JW, Verjee Y, Hill EW, Rege
483–493. JEO (2002). African pastoralism: genetic imprints of origins
Dray S, Dufour AB (2007). The ade4 package: implementing the and migrations. Science 296: 336–339.
duality diagram for ecologists. J Stat Softw 22: 1–20. Harville DA (1997). Matrix Algebra From a Statistician’s Perspec-
Dray S, Chessel D, Thioulouse J (2003a). Co-inertia analysis and tive. Springer: New York.
the linking of ecological data tables. Ecology 84: 3078–3089. Hotelling H (1936). Relations between two sets of variates.
Dray S, Chessel D, Thioulouse J (2003b). Procrustean co-inertia Biometrika 28: 321–327.
analysis for the linking of multivariate datasets. Ecoscience 10: Jambu M (1991). Exploratory and Multivariate Data Analysis.
110–119. Academic Press Inc.: Orlando, Florida.
Dray S, Dufour AB, Chessel D (2007). The ade4 package—II: Jarraud S, Mougel C, Thioulouse J, Lina G, Meugnier H, Forey F
Two-table and K-table methods. R News 7: 47–54. et al. (2002). Relationships between Staphylococcus aureus
Dray S (2008). On the number of principal components: A test genetic background, virulence factors, agr groups (alleles),
of dimensionality based on measurements of similarity and human disease. Infect Immun 70: 631–641.
between matrices. Comput stat data anal 52: 2228–2237. Johnson FM, Schaffer HE, Gillaspy JE, Rockwood ES (1969).
Dray S, Saı̈d S, Debias F, Chessel D (2008). Spatial ordination of Isozyme genotype-environment relationships in natural
vegetation data using a generalization of Wartenberg’s populations of the harvester ant, Pogonomyrmex barbatus,
multivariate spatial correlation. J Veg Sci 19: 45–56. from Texas. Biochem Genet 3: 429–450.
Dupanloup I, Schneider S, Excoffier L (2002). A simulated Johnson FM, Schaffer HE (1973). Isozyme variability in species
annealing approach to define the genetic structure of of the genus drosophila. VII. Genotype-environment rela-
populations. Mol Ecol 11: 2571–2581. tionships in populations of D. melanogaster from the eastern
Edwards AWF (2003). Human genetic diversity: Lewontin’s United States. Biochem Genet 10: 149–163.
fallacy. BioEssays 25: 798–801. Jombart T (2008). adegenet: a R package for the multivariate
Escoufier Y (1987). The duality diagramm: a means of better analysis of genetic markers. Bioinformatics 24: 1403–1405.
practical applications. In: Legendre P, Legendre L (eds). Jombart T, Devillard S, Dufour AB, Pontier D (2008). Revealing
Development in Numerical Ecology. NATO advanced Institute, cryptic spatial patterns in genetic variability by a new
Serie G. Springer Verlag, Berlin. pp 139–156. multivariate method. Heredity 101: 92–103.

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

340
Krzanowski WJ, Marriott FHC (1995). Multivariate Analysis. Part Patterson N, Price AL, Reich D (2006). Population structure and
2: Classification, Covariance Structures and Repeated Measure- eigenanalysis. PLoS genet 2: 2074–2093.
ments. Halsted Press, John Wiley & Sons: Edward Arnold, Pavoine S, Dufour AB, Chessel D (2004). From dissimilarities
London. among species to dissimilarities among communities:
Kölliker R, Bassin S, Schneider D, Widmer F, Fuhrer J (2008). a double principal coordinate analysis. J Theor Biol 228:
Elevated ozone affects the genetic composition of Plantago 523–537.
lanceolata L. Populations. Environ Pollut 152: 380–386. Pavoine S, Bailly X (2007). New analysis for consistency among
Lachenbruch PA, Goldstein M (1979). Discriminant analysis. markers in the study of genetic diversity: development and
Biometrics 35: 69–85. application to the description of bacterial diversity. BMC
Laloë D, Jombart T, Dufour AB, Moazami-Goudarzi K (2007). Evolut Biol 7: 156.
Consensus genetic structuring and typological value of Pearson K (1901). On lines and planes of closest fit to systems of
markers using multiple co-inertia analysis. Genet Sel Evol points in space. Philos Mag 2: 559–572.
39: 545–567. Perrière G, Thioulouse J (2002). Use and misuse of correspon-
Lebart L, Morineau A, Piron M (2004). Statistique Exploratoire dence analysis in codon usage studies. Nucleic Acids Res 30:
Multidimensionnelle. DUNOD: Paris. 4548–4555.
Legendre P, Legendre L (1998). Numerical Ecology. Elsevier Powell JR, Taylor CE (1978). Are human races ‘substantially’
Science B.V.: Amsterdam. different genetically? Am Nat 112: 1139–1142.
Legendre P, Anderson DJ (1999). Distance-based redundancy Preziosi RF, Fairbairn DJ (1992). Genetic population structure
analysis: testing multispecies responses in multifactorial and levels of gene flow in the stream dwelling waterstrider
ecological experiments. Ecol Monogr 69: 1–24. Aquarius ( ¼ Gerris) remigis (Emiptera: Geridae). Evolution 46:
Lessa EP (1990). Multidimensional analysis of geographic 430–444.
genetic structure. Syst Zool 39: 242–252. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA,
Lewontin RC (1978). Single-locus and multiple-locus Reich D (2006). Principal components analysis corrects for
measures of genetic distance between groups. Am Nat 112: stratification in genome-wide association studies. Nat Genet
1138–1139. 38: 904–909.
Li MH, Zhao SH, Bian C, Wang HS, Wei H, Liu B et al. (2002). Pritchard JK, Stephens M, Donnelly P (2000). Inference of
Genetic relationships among twelve chinese indigenous goat population structure using multilocus genotype data.
populations based on microsatellite analysis. Genet Sel Evol Genetics 155: 945–959.
34: 729–744. R Development Core Team (2008). R: A Language and Environ-
MacHugh DE, Shriver MD, Loftus RT, Cunningham P, ment for Statistical Computing. R Foundation for Statistical
Bradley DG (1997). Microsatellite DNA variation and the Computing: Vienna, Austria ISBN 3-900051-07-0.
evolution, domestication and phylogeography of taurine http://www.R-project.org.
and zebu cattle (Bos taurus and Bos indicus). Genetics 146: Rao CR (1964). The use and interpretation of principal
1071–1086. component analysis in applied research. Sankhya, A 26:
MacHugh DE, Loftus RT, Cunningham P, Bradley DG 329–359.
(1998). Genetic structure of seven European cattle Reyment RA (2005). The statistical analysis of multivariate
breeds assessed using 20 microsatellite markers. Anim Genet serological frequency data. Bull Math Biol 67: 1303–1313.
29: 333–340. Rohlf FJ (1963). Classification of Aedes by numerical
Matsuoka Y, Vigouroux Y, Goodman MM, Jesus Sanchez G, taxonomic methods (diptera:culicidae). Ann Entomol Soc Am
Buckler E, Doebley J (2002). A single domestication for maize 56: 798–804.
shown by multilocus microsatellite genotyping. Proc Natl Romano V, Calı́ F, Ragalmuto A, D’Anna RP, Flugy A, De Leo G
Acad Sci USA 99: 6080–6084. et al. (2003). Autosomal microsatellite and mtDNA genetic
McKechnie SW, Ehrlich PR, White RR (1975). Population analysis in Sicily (Italy). Ann Hum Genet 67: 42–53.
genetics of euphydryas butterflies. I. genetic variation and Sagnard F, Barberot C, Fady B (2002). Structure of
the neutrality hypothesis. Genetics 81: 571–594. genetic diversity in Abies alba Mill. from southwestern
McRae BH, Beier P, Dewald LE, Huynh LY, Keim P (2005). Alps: multivariate analysis of adaptive and non-adaptative
Habitat barriers limit gene flow and illuminate historical traits for conservation in France. For Ecol Manage 157:
events in a wide-ranging carnivore, the American puma. 175–189.
Mol Ecol 14: 1965–1977. Sanchez-Mazas A, Langaney A (1988). Common genetic pools
Menozzi P, Piazza A, Cavalli-Sforza LL (1978). Synthetic between human populations. Hum Genet 78: 161–166.
maps of human gene frequencies in Europeans. Science 201: Schaffer HE, Johnson FM (1974). Isozyme allelic frequencies related
786–792. to selection and gene-flow hypotheses. Genetics 77: 163–168.
Mitton JB (1978). Measurement of differentiation: reply to Seal HL (1966). Multivariate Statistical Analysis for Biologists.
Lewontin, Powell and Taylor. Am Nat 112: 1142–1144. Methuen and co.: London.
Moazami-Goudarzi K, Laloë D, Furet JP, Grosclaude F (1997). Seber GAF (1977). Linear Regression Analysis. John Wiley & Sons:
Analysis of genetic relationships between 10 cattle breeds New York.
with 17 microsatellites. Anim Genet 28: 338–345. She JX, Autem M, Kotulas G, Pasteur N, Bonhomme F (1987).
Moazami-Goudarzi K, Laloë D (2002). Is a multivariate Multivariate analysis of genetic exchanges between Solea
consensus representation of genetic relationships among aegyptiaca and Solea senegalensis (Teleosts, Soleidae). Biol J
populations always meaningful? Genetics 162: 473–484. Linnean Soc 32: 357–371.
Mulley JC, James JW, Barker JSF (1979). Allozyme genotype- Smouse PE, Spielman RS, Park MH (1982). Multiple-locus
environment relationships in natural populations of allocation of individuals to groups as a function of the
Drosophila buzzatii. Biochem Genet 17: 105–126. genetic variation within and differences among human
Pariset L, Savarese MC, Cappuccio I, Valentini A (2003). Use of populations. Am Nat 119: 445–463.
microsatellites for genetic variation and inbreeding analysis Soshnikov A, Fyodorov YV (2005). On the largest singular
in Sarda sheep flocks of central Italy. J Anim Breed Genet 120: values of random matrices with independent Cauchy entries.
425–432. J Math Phys 46: 033302.
Parisod C, Christin PA (2008). Genome-wide association to Takeuchi K, Yanai H, Mukherjee BN (1982). The foundations of
fine-scale ecological heterogeneity within a continuous multivariate analysis: a unified approach by means of
population of Biscutella laevigata (brassicaceae). New Phytol projection onto linear subspaces. Wiley Eastern Limited:
178: 436–447. New-Delhi.

Heredity
 Multivariate analysis of genetic markers
T Jombart et al

341
Taylor CE, Mitton JB (1974). Multivariate analysis of genetic Weir BS (1996). Genetic Data Analysis II. Sinauer Associates:
variation. Genetics 76: 575–585. Sunderland, Massachussetts.
Ter Braak CJF (1986). Canonical correspondence analysis : a Williams BK, Titus K (1988). Assessment of sampling stability in
new eigenvector technique for multivariate direct gradient ecological applications of discriminant analysis. Ecology 69:
analysis. Ecology 67: 1167–1179. 1275–1285.
van Pijlen IA, Amos B, Burke T (1995). Patterns of genetic Xuebin Q, Jianlin H, Chekarova I, Badamdorj D, Rege JEO,
variability at individual minisatellite loci in minke whale Hanotte O (2005). Genetic diversity and differentiation of
Balaenoptera acutorostrata populations from three different Mongolian and Russian yak populations. J Anim Breed Genet
oceans. Mol Biol Evol 12: 459–472. 122: 117–126.
Warnes GR (2003). The genetics package. R News 3: 9–13. Zhivotovsky LA, Rosenberg NA, Feldman MW (2003). Features
Wartenberg DE (1985). Canonical trend surface analysis: of evolution and expansion of modern humans, inferred
a method for describing geographic patterns. Syst Zool 34: from genomwide microsatellite markers. Am J Hum Genet 72:
259–279. 1171–1186.

Heredity