## Abstract

With the advent of modern imaging and measurement technology, complex phenotypes are increasingly represented by large numbers of measurements, which may not bear biological meaning one by one. For such multivariate phenotypes, studying the pairwise associations between all measurements and all alleles is highly inefficient and prevents insight into the genetic pattern underlying the observed phenotypes. We present a new method for identifying patterns of allelic variation (genetic latent variables) that are maximally associated—in terms of effect size—with patterns of phenotypic variation (phenotypic latent variables). This multivariate genotype–phenotype mapping (MGP) separates phenotypic features under strong genetic control from less genetically determined features and thus permits an analysis of the multivariate structure of genotype–phenotype association, including its dimensionality and the clustering of genetic and phenotypic variables within this association. Different variants of MGP maximize different measures of genotype–phenotype association: genetic effect, genetic variance, or heritability. In an application to a mouse sample, scored for 353 SNPs and 11 phenotypic traits, the first dimension of genetic and phenotypic latent variables accounted for >70% of genetic variation present in all 11 measurements; 43% of variation in this phenotypic pattern was explained by the corresponding genetic latent variable. The first three dimensions together sufficed to account for almost 90% of genetic variation in the measurements and for all the interpretable genotype–phenotype association. Each dimension can be tested as a whole against the hypothesis of no association, thereby reducing the number of statistical tests from 7766 to 3—the maximal number of meaningful independent tests. Important alleles can be selected based on their effect size (additive or nonadditive effect on the phenotypic latent variable). This low dimensionality of the genotype–phenotype map has important consequences for gene identification and may shed light on the evolvability of organisms.

STUDIES of genotype–phenotype association are central to several branches of contemporary biology and biomedicine, but they suffer from serious conceptual and statistical problems. Most of these studies consist of a vast number of pairwise comparisons between single genetic loci and single phenotypic variables, typically leading—among other reasons—to very low fractions of phenotypic variance explained by genetic effects [“missing heritability” (Manolio *et al.* 2009; Eichler *et al.* 2010)]. *Post hoc* corrections for multiple testing can lead to a dramatic loss of statistical power and in fact violate standard rules of statistical inference. Biologically more important, most phenotypes are not determined by single alleles, but by the joint effects, both additive and nonadditive, of a number of alleles at multiple loci. With the advent of modern imaging and measurement technology, complex phenotypes, such as the vertebrate brain or cranium, often are represented by large numbers of variables. This further complicates the study of genotype–phenotype association by tremendously increasing the number of pairwise comparisons between genetic loci and phenotypic variables, which may not be meaningful traits *per se* [for instance, in geometric morphometrics, voxel-based image analysis, and many behavioral studies (Bookstein 1991; Ashburner and Friston 2000; Mitteroecker and Gunz 2009; Houle *et al.* 2010)]. The genotype–phenotype associations we actually seek are between certain allele combinations from multiple loci and certain combinations of phenotypic variables that bear biological interpretation. The number of such pairs of “latent” allele combinations and phenotypes that underlie the observed genotype–phenotype association depends on the genetic-developmental system under study, but typically is less than the number of assessed loci and phenotypic variables (Hallgrimsson and Lieberman 2008; Martinez-Abadias *et al.* 2012).

Several methods have been suggested for such a multivariate mapping, including multiple and multivariate regression (Haley and Knott 1992; Jansen 1993; Hackett *et al.* 2001; de Los Campos *et al.* 2013), principal component regression (Wang and Abbott 2008), low-rank regression models (Zhu *et al.* 2014), partial least-squares regression (Bjørnstad *et al.* 2004; Bowman 2013), and canonical correlation analysis (Leamy *et al.* 1999; Ferreira and Purcell 2009). We present a multivariate analytic strategy—which we term *multivariate genotype–phenotype mapping* (MGP)—that embraces and relates all of these methods and that circumvents several of the problems resulting from pairwise univariate mapping and from the multivariate analysis of the loci separately from the phenotypes. Our approach does not primarily aim for the detection and location of single loci segregating with a given phenotypic trait. Instead, we present an approach that identifies patterns of allelic variation that are maximally associated—in terms of effect size—with patterns of phenotypic variation. In this way, we gain insight into the multivariate structure of genotype–phenotype association, including its dimensionality and the clustering of genetic and phenotypic variables within this association—the genetic-developmental properties determining the evolvability of organisms (Wagner and Altenberg 1996; Hendrikse *et al.* 2007; Mitteroecker 2009; Pavlicev and Hansen 2011).

## The Principle of Multivariate Genotype–Phenotype Mapping

Let there be *p* genetic loci and *q* phenotypic measurements scored for *n* specimens. Instead of assessing each of the pairwise genotype–phenotype associations, we seek a genetic effect—composed of the additive and nonadditive effects of multiple alleles—onto a phenotypic trait that is a composite of multiple measured phenotypic variables. As these genetic and phenotypic features are not directly measured, but perhaps present in the data, we refer to them as genetic and phenotypic “latent variables,” and (Figure 1). The molecular, physiological, and developmental processes that underlie the genotype–phenotype relationship and that constitute the latent variables likely are complex nonlinear processes. In a first approximation, however, we consider the latent variables as linear combinations of the alleles and phenotypic measurements, respectively.

How do we identify these latent variables? This problem can be viewed in a dual way. First, it can be assumed that the effect of a genetic latent variable on a phenotypic latent variable is stronger than the effect of any single locus on any single phenotypic variable. We may thus seek genetic and phenotypic latent variables (linear combinations) with *maximal genotype–phenotype association*. In addition, there might be further pairs or “dimensions” of genetic and phenotypic latent variables with maximal associations, mutually independent across the dimensions, that together account for the observed genotype–phenotype association. The classic measures of genotype–phenotype association in the quantitative genetic literature are (i) genetic effect (average, additive, and dominance effect); (ii) genetic variance (the phenotypic variance accounted for by genetic effects); and (iii) heritability, the ratio of genetic to total phenotypic variance. Accordingly, we may compute latent variables that maximize one of these measures of genotype–phenotype association, depending on the scientific question and the information content in the data.

A second, equivalent, way to view the problem is that of a search for simple patterns underlying the observed pairwise genotype–phenotype associations. Technically, we seek low-rank (*i.e.*, “simple”) matrices that approximate the matrix of pairwise associations. A powerful standard technique in multivariate statistics for this purpose is singular value decomposition (SVD). For a matrix of pairwise genotype–phenotype associations, SVD finds pairs of singular vectors (one of length *p*, one of length *q*) of which the outer product best approximates the matrix in a least-squares sense. The two singular vectors can be interpreted as the genetic and phenotypic effects of the corresponding latent variables (the coefficients in Figure 1) that best approximate the observed genotype–phenotype associations. When connecting the two views, effect maximization and pattern search, the question arises: For which kind of matrices do the singular vectors (as simple patterns) lead to latent variables that maximize the above measures of genotype–phenotype association?

In *Materials and Methods* we demonstrate how to identify these latent genetic and phenotypic variables by an SVD of the appropriate association matrix. A more detailed derivation, including proofs, is given in *Appendix A*. In an application to a classic mouse sample, scored for 353 genetic markers and 11 phenotypic variables, we demonstrate the effectiveness of multivariate genotype–phenotype mapping. We show that a single dimension of latent variables already suffices to account for >70% of genetic variance present in the 11 traits. The first three dimensions together account for almost 90% of genetic variation and capture all the interpretable genotype–phenotype association in the data. Each dimension can be tested as a whole against the hypothesis of no association by a permutation approach, which reduces the number of significance tests from 7766 to 3. We discuss the consequences of this low dimensionality of the genotype–phenotype map for gene identification and for understanding the evolvability of organisms.

## Materials and Methods

### Maximizing genotype–phenotype association

Let each allele at each locus be represented by a vector that contains the additive genotype scores (0, 1, or 2) for all *n* subjects. Additional vectors may represent interactions of alleles at one locus (dominance scores, 0 or 1). See below for the implementation of epistasis, the interaction of alleles at different loci. We seek the combined effect of the alleles on a phenotype composed of the measured variables . For notational convenience, let the columns of **X** and **Y** be mean centered. Let the effects of the alleles on the phenotype be denoted by the vector and the weightings of the measured variables that determine the phenotype by the vector . Then the genetic and phenotypic latent variables are given by the linear combinations and (Figure 1). Since, by definition, the latent variables have a stronger genotype–phenotype association than any single variable, the coefficient vectors are chosen to maximize the association between the corresponding latent variables: (i) genetic effect, (ii) genetic variance, or (iii) heritability (see *Appendix A*, *Univariate genotype–phenotype* *association*). In addition to this pair of latent variables, there are further pairs of latent variables with effects and independent of the previous ones, that together account for the observed genotype–phenotype association.

For any real matrix, SVD yields a first pair of real singular vectors both of unit length, and a real singular value The “left” singular vector is of dimension and the “right” vector is of dimension The outer product scaled by is the rank-1 matrix that best approximates the matrix in a least-squares sense. There is also a second pair of singular vectors orthogonal to the first singular vectors, which are associated with a second singular value Together, the two pairs of singular vectors give the best rank-2 approximation of the matrix, and so forth for further dimensions. SVD yields optimal low-rank approximations of the original matrix by *maximizing* the singular values, that is, the contribution of the corresponding rank-1 matrix to the matrix approximation. The summed squared singular values equal the summed squared elements of the approximated matrix. The number of relevant dimensions (the rank of the approximation) can thus be determined by the squared singular values, expressed as a fraction of the total squared singular values.

For a matrix of pairwise genotype–phenotype associations, the singular vectors may serve as weightings for the genetic and phenotypic variables to compute the latent variables and The question is, For which kind of association matrices are the singular vectors the vectors maximizing (i) genetic effect, (ii) genetic variance, and (iii) heritability, respectively?

#### Genetic effect (i):

The *additive genetic effect* of one allele substitution is half the difference between the homozygote mean phenotypes, and the *dominance effect* is the deviation of the heterozygote mean phenotype from the midpoint of the homozygote mean phenotypes. By contrast, the *average effect* of an allele substitution is the average difference between offspring that get this allele and random offspring (Falconer and Mackay 1996; Roff 1997). In a sample of measured individuals, the average effect can be estimated by the regression slope of the phenotype on the additive genotype scores, whereas additive and dominance effects can be estimated by the two multiple-regression coefficients of the phenotype on both the additive and dominance genotype scores (see *Appendix A*, *Univariate genotype–phenotype* *association*).

For multivariate genotypes and phenotypes, maximizing the effect of the genetic latent variable on the phenotypic latent variable translates into finding vectors **a** and **b** that maximize the slope (1)of the regression of the phenotype on the allele combination Under the constraint this regression slope is maximized by the first pair of singular vectors of the matrix of regression coefficients (2)For detailed derivations and proofs see *Appendix A*, *Multivariate genotype–phenotype association* and *Maximizing genotype–phenotype association via SVD*. See *Appendix* *B*, *Computational properties* for a discussion of the computational difficulties arising from the matrix inversion if

If the genetic variables **X** comprise the *p* additive genotype scores only, maximizing the regression slope (1) maximizes the average effect, and the linear combination can be interpreted as breeding values (sum of average effects). The *p* elements of **a** are the partial average effects of the corresponding alleles on the phenotype (average effects conditional on the other alleles); they are proportional to the regression coefficients of on **X**. The maximal slope (maximal average effect) associated with the linear combinations is given by the singular value

If both additive and dominance scores are included in **X**, the sum of both additive and dominance effects is maximized and the linear combination can be interpreted as genotypic values. The elements of **a**, which is now of dimension would then correspond to the partial additive and dominance effects on the phenotype and the singular value is the sum of additive and dominance effects.

#### Genetic variance (ii):

For a population in Hardy–Weinberg equilibrium, the additive genetic variance can be expressed as the product of the variance of the additive genotype scores and the squared average effect (*cf*. Equation A3 in *Appendix A*). Maximizing genetic variance for multivariate genotypes and phenotypes thus is achieved by maximizing over the vectors **a** and **b**. The latent variables with maximal genetic variance are given by the vectors and where and are the first left and right singular vectors of the matrix (3)If the genetic variables comprise the additive scores only, this approach maximizes the additive genetic variance, given by whereas if both additive and dominance scores are included, the total genetic variance (additive plus dominance variance) is maximized. This approach is computationally equivalent to reduced-rank regression (Izenman 1975; Aldrin 2000).

#### Heritability (iii):

Heritability can be expressed as the squared correlation coefficient of the phenotype and the genotype scores (*Appendix A*, *Univariate genotype–phenotype association*). Hence, maximizing heritability in the multivariate context amounts to maximizing the squared correlation which is achieved by the vectors and where and are the left and right singular vectors of the matrix (4)Including only the additive genotype scores in **X** maximizes the narrow-sense heritability whereas including both additive and dominance scores maximizes the broad-sense heritability resulting from additive and dominance variance. These maximal heritabilities equal the squared singular value This approach is equivalent to canonical correlation analysis (*e.g.*, Mardia *et al.* 1979).

#### Covariance (iv):

Bjørnstad *et al.* (2004) and Mehmood *et al.* (2011) applied partial least-squares analysis (PLS) to identify genetic and phenotypic latent variables. PLS maximizes the covariance between the two linear combinations The unit vectors maximizing this covariance can be computed as the first pair of singular vectors of the cross-covariance matrix The maximal covariance is given by The covariance between a phenotypic variable and a genetic variable has no correspondence in the classic genetic framework. However, this approach has convenient computational properties as it requires no matrix inverse (see *Appendix B*, *Computational properties*). The scaled genetic coefficients are equal to the covariances between the corresponding locus and the phenotypic latent variable, without conditioning on the other loci as in approaches i–iii. See *Appendix A*, *Partial* *least-squares analysis* for more details.

### Properties of the four approaches

For each of the four approaches (Table 1), further dimensions (pairs of latent variables ) can be extracted by the subsequent pairs of singular vectors of the corresponding association matrix. The second dimension consists of a new allele combination (genetic latent variable) and a new phenotype (phenotypic latent variable) that are independent of the ones from the first dimension and have the second largest association, and similarly for further dimensions. The maximal number of dimensions is However, the notion of “independence” differs among the four approaches. In approaches i and iv—the maximizations of genetic effect and covariance—the genetic coefficient vectors are orthogonal ( for and 1 for ), whereas in approaches ii and iii—the maximizations of genetic variance and heritability—the genetic latent variables are uncorrelated ( for and 1 for ). The phenotypic effects are orthogonal in approaches i, ii, and iv, whereas the phenotypic latent variables are uncorrelated in approach iii.

In the first three approaches, the corresponding association is maximized conditional on all other linear combinations of alleles, including the other latent variables. The matrix of regression coefficients of the phenotypic latent variables on the genetic latent variables thus is diagonal, where the matrix **A** contains the vectors and **B** the vectors In words, the prediction of the *i*th phenotypic latent variable by the *i*th genetic latent variable is not improved by adding any other latent variable or any linear combination of them as predictor. In this sense, the effect of the allele combination on the phenotype is “independent” of that of any other allele combination Furthermore, in approaches ii, iii, and iv the genetic latent variables are correlated only with the corresponding phenotypic latent variable but not with any of the other latent variables: for (see *Appendix A*, *Maximizing genotype–phenotype association via SVD* for proofs).

The constraints on the mutual independence (orthogonality or uncorrelatedness) of the genetic coefficient vectors and the phenotypic coefficient vectors are unlikely to reflect biological relationships and can complicate the reification of multiple latent variables as biological factors (see, *e.g.*, Cheverud 2007; Bookstein 2014). Supplemental Material, Figure S1 illustrates these properties by applications to simple simulated data sets. Unless only a single dominant dimension (pair of latent variables) is present in the data or all the important latent variables can clearly be interpreted, multiple dimensions should be interpreted jointly, as spanning a genetic subspace that is maximally associated with a phenotypic subspace. The summed squared singular values of these dimensions indicate their joint genotype–phenotype association.

Maximizing genetic effect (phenotypic change per unit genetic change) in approach i imposes a constraint on the norms of **a** and **b**; that is, it implies a concept of “total” genetic and phenotypic effect. The vectors are computed to have a norm (summed squared elements) of 1, which is common in statistics but not an obvious choice in genetics (see *Appendix B*, *Geometric properties* for further discussion). By contrast, maximizing genetic variance and heritability in approaches ii and iii requires no constraint on the norm of **a**; the explained phenotypic variance is unaffected by the norm of **a**. (The SVD scales **a** so that the latent variable has unit variance, but this choice has no effect on the maximal genetic variance.) In approaches ii and iii, the singular values and the right singular vectors are even unchanged under all linear transformations of the genetic variables **X** (see *Appendix B*, *Geometric properties* for a proof and Figure S1 for a demonstration). This implies that the maximal genetic variances (singular values) in approach ii as well as the phenotypes that show these maximal variances (right singular vectors) do not depend on the variances and covariances of the genetic variables (as they are scale dependent), that is, on genetic variance and linkage disequilibrium. They depend on the *relative* differences between the genotype scores only. The same applies to approach iii, but the heritabilities are even invariant to affine transformation of the phenotypic variables. To understand these properties intuitively, consider that the regression slope (genetic effect) depends on the variance of the predictor variable (the genotype scores), but the explained variance (genetic variance, heritability) does not. The genetic coefficients (left singular vectors) are affected by linear transformations of **X** because they are the multiple-regression coefficients of the phenotypic latent variable on the loci.

Approaches i, ii, and iv require meaningful phenotypic variances and covariances that are commensurate across the variables (because **b** is constrained to a norm of 1). This can be the case in morphometrics and chemometrics, but rarely in the behavioral sciences and psychometrics (see also Mitteroecker and Huttegger 2009; Huttegger and Mitteroecker 2011). If correlations, but not variances, are interpretable, for example when the variables have incommensurate scales, approaches i, ii, and iv may be applied after standardizing the variables separately to unit variance or unit mean (Hansen and Houle 2008). However, linear combinations of variance-standardized variables may be difficult to interpret; multivariate genotype–phenotype mapping is most powerful in contexts such as geometric morphometrics and image analysis, where the variables do not require standardization and interpretations are always in terms of linear combinations of measurements. When neither variances nor correlations are interpretable for the phenotypes, only approach iii—the maximization of heritability—can be applied.

### Epistasis

The interactions of alleles at different loci are not considered explicitly in the above descriptions. Epistatic effects contribute to the estimates as far as they contribute to the average or additive allelic effects (Cheverud and Routman 1995, 1996). However, it is possible to explicitly include epistatic effects by appending a design matrix to **X** that represents pairwise or higher-order allele interactions. When adding such interaction terms, approach i maximizes the sum of additive, dominance, and epistatic effects; approach ii maximizes genetic variance resulting from all three kinds of genetic effects (additive plus interaction variance); and approach iii maximizes the corresponding heritability.

As noted earlier, the allele combinations have no mutual interactions in their effects on the corresponding phenotypic latent variable and in this sense show no epistasis. But this property should not be overinterpreted; it is primarily a convenient mathematical property, which holds true for all data sets.

### Statistical significance

The presented approaches are exploratory techniques for identifying multivariate patterns underlying the observed genotype–phenotype relationships. As long as no specific hypotheses about the effect of particular alleles have been formulated, no allele-specific hypothesis tests are appropriate. But one can test each pattern (pair of genetic and phenotypic latent variables) as a whole against the hypothesis of no genotype–phenotype association by a permutation test (*e.g.*, Churchill and Doerge 1994; Good 2000): Randomly reassign the phenotype vector of each individual to another individual in the data set and recompute the singular value. The *P*-value is given by the fraction of such permutations that yield an equal or larger singular value than the original one. A structured pedigree in an experimental cross may require a restricted or hierarchical permutation approach. Confidence intervals of the genetic and phenotypic coefficients can be estimated by a bootstrap approach (*e.g.*, Visscher *et al.* 1996). These intervals can be used to compare the computational stability of coefficients across alleles or across dimensions. For the reasons mentioned earlier, they should not be used for inferences about allele-specific statistical significance.

### Data and analysis

We demonstrate multivariate genotype–phenotype mapping by an application to a classic mouse sample, comprising 1581 specimens of the F_{2} and F_{3} generations of an intercross of inbred LG/J and SM/J mice, obtained from The Jackson Laboratory (for details see Cheverud *et al.* 1996; Norgard *et al.* 2008). Each individual was genotyped at 384 polymorphic SNPs, of which we used the 353 SNPs on the 19 autosomal chromosomes. Every locus was represented by one variable for the additive genotype scores and one variable for the dominance scores. Total body weight (WTN) and the weights of the reproductive fat pad (FP), the heart (HT), the kidney (KD), the spleen (SP), and the liver (LV), as well as tail length (TL), were recorded for each individual. In addition, the lengths of the right set of long bones (humerus, ulna, femur, and tibia) were measured. Measurements were corrected for the effects of sex, litter size, and age at necropsy and standardized to unit variance. All work on mice was performed using a protocol approved by the Institutional Animal Care and Use Committee of the Washington University School of Medicine (St. Louis). These data may not be representative of modern high-throughput standards, but they allowed us to test our approach on a well-explored sample with easily interpretable variables.

We analyzed these data by approach ii—the maximization of total genetic variance—because genetic variance is the focus of most breeding studies and evolutionary models. The computation of the inverse square root matrix was based on a generalized inverse, using the first 80 principal components of the genetic variables, which accounted for 76% of total variation (see *Appendix B*, *Computational* *properties*). To check for overfitting the data (see *Discussion* and *Appendix B*, *Computational* *properties*), we used a leave-one-out cross-validation: Each individual’s phenotype was predicted by the genetic and phenotypic coefficient vectors that were computed using only the other individuals. The explained variance was then computed from these predictions.

### Data availability

All data used in the analysis are available at the DRYAD repository: DOI:10.5061/dryad.fc55k. All computations and visualizations were performed in Mathematica 9.0 (Wolfram Research Inc.).

## Results

Figure 2A shows the genetic variances for the 11 dimensions of latent variables in the style of a scree plot. The first dimension strikingly dominates the genotype–phenotype map by accounting for 72% of the total genetic variance summed over all 11 phenotypic variables (Table 2). The corresponding phenotypic latent variable represented overall limb length (high positive coefficients for all long bones; Figure 3A). More than 43% of the variance in this phenotypic latent variable was accounted for by the corresponding genetic latent variable (Table 2). The additive and dominance effects of the loci are represented by the genetic coefficients shown in Figure 4. They indicate strong additive effects on limb length mainly on chromosomes 1, 2, 3, 6, 8, 9, 13, and 17.

The second phenotypic latent variable reflected body and organ weight (Figure 3B), with additive and dominance effects mainly on chromosomes 2, 3, 6, 9, 10, 11, 12, 13, and 19 (Figure S2). Compared to dimension 1, the explained variance of body/organ weight was relatively small (7% of phenotypic variance). This pattern differs from the third phenotypic latent variable, which contrasted distal *vs.* proximal long bone length. This trait varied little in the studied population but was under stronger genetic control than body weight (25% explained phenotypic variance). The strongest additive genetic effects were located on chromosomes 1 and 11 (Figure S3). These three dimensions together accounted for 87% of the genetic variance present in the 11 variables of this sample. The subsequent dimensions accounted for small portions of genetic variance only (<3%) and had no obvious interpretation.

Figure 5 shows a plot of the three scaled phenotypic coefficient vectors, constituting a phenotype space in which the phenotypic variables cluster according to their genetic structure. For approach ii applied here, the squared length of the vectors approximates the genetic variance of the corresponding phenotypic variable, and the cosine of the angle between the vectors approximates their genetic correlation. The long bones clustered together to the exclusion of the weight measurements, indicating a shared genetic basis. Humerus and femur as well as tibia und ulna showed particularly strong genetic correlations. Tail length was more closely correlated with the weight measurements than with the long bone lengths.

In a second analysis, we performed a separate and more detailed study of chromosome 6, which showed strong genetic effects. As genetic predictors we used the additive and dominance genotype scores for each of the 22 screened SNPs as well as for 89 loci imputed every 1 cM. This allows for the implementation of interval mapping (Lander and Botstein 1989) in multivariate genotype–phenotype mapping. These 222 genetic variables and the 11 phenotypic variables were analyzed again with approach ii, the maximization of genetic variance, resulting in two interpretable dimensions (Figure 2B and Figure 6). Bootstrap confidence intervals are shown for both the genetic and phenotypic effects. The computation of the matrix inverse was based on the first eight principal components of the genetic variables (89% of variance). The first dimension again represented limb length and accounted for even 85% of total genetic variance within the 11 phenotypic variables (). The additive genetic effects had two distinguished peaks, one at ∼60–70 Mb and one at ∼140 Mb. Both peaks were associated with small dominance effects. These results correspond well to the two loci affecting long bone length identified by Norgard *et al.* (2008); they were estimated at 85 and 144 Mb on chromosome 6 by applying traditional methods to the same data. The second dimension represented body and organ weight and accounted for 8% of total genetic variance (). Additive and dominance effects had three peaks, one close to the centromere, one at ∼60–70 Mb, and one at the end of the chromosome. For the latter region, additive and dominance effects were of opposite sign. These three locations are in accordance with the loci identified by Vaughn *et al.* (1999) and Fawcett *et al.* (2008) for body and organ weight in the F_{2} and F_{3} generations of the same cross. The confidence intervals of both genetic and phenotypic coefficients for dimension 2 were considerably wider than that for dimension 1. Note that the confidence intervals should not be used for statistical inference in this exploratory context, only for the comparison of computational stability within the sample.

Since the F_{3} population consisted of 200 sets of full sibs, we accounted for the genetic similarity between individuals in a separate analysis. We constructed an matrix that represents genetic similarity between individuals (we tested both the expected similarity based on relatedness and the actual similarity based on the SNP data) and implemented this matrix in the estimation by a generalized least-squares approach (see *Appendix B*, *Generalized least squares*). This had only a limited effect on the results and we thus presented just the ordinary least-squares solutions here. We also repeated the analyses with different numbers of principal components (PCs) used for the matrix inversion. The genetic coefficients were stable against small changes of the numbers of PCs; the phenotypic coefficients and the shape of the scree plot were stable even over a very wide range of PCs. Taking the cube root of the weight measurements before variance standardization had basically no effect on the results. We checked whether outliers could drive some of the results, but found no evidence in scatter plots of genetic *vs.* phenotypic latent variables. We also applied the other three maximization approaches to the data, which basically resulted in the same first two pairs of phenotypic and genetic latent variables with similar explained variances (see *Appendix C* for a brief presentation of these results).

## Discussion

Many traits are affected by numerous alleles with small or intermediate effects, which are difficult to detect by mapping each locus separately. Loci selected by separately computed *P*-values often account for low fractions of phenotypic variance and provide an incomplete picture of the genotype–phenotype map (Manolio *et al.* 2009; Eichler *et al.* 2010). “Whole-genome prediction” methods, which are based on all scored loci, have proved more effective in explaining phenotypic variation of a trait (Yang *et al.* 2010; de Los Campos *et al.* 2013). However, many complex phenotypes cannot be adequately represented by a single variable but require multiple measurements. For such multivariate traits, we presented a “whole-genome, whole-phenotype prediction” method that identifies the genetically determined traits and their associated allelic effects in a single step. This avoids an inefficient decomposition (such as principal component analysis) of the phenotypic variables separately from the genetic variables (Cheverud 2007); instead our method provides a decomposition of the genotype–phenotype map itself.

Multivariate genotype–phenotype mapping separates phenotypic features under strong genetic determination from features with less genetic control, whereas traditional mapping of complex traits typically lumps different phenotypic features with different heritabilities. For example, we found that overall limb length in our mouse sample is both highly variable and highly heritable (43% explained phenotypic variance); a second feature, distal *vs.* proximal limb bone length, is much less variable but still shows an explained variance of 25%. All other aspects of long bone variation, *e.g.*, forelimb *vs.* hindlimb length, show very little genetic variation. Mapping each of the four limb bones separately thus leads to estimates of explained variance that are averages across all these features, some of which have high heritability and some of which have basically none. It thus misses the actual signal: the traits (latent variables) under strong genetic control. Multivariate genotype–phenotype mapping can tell one *where* to look for the (missing) heritability in the phenotype and shows the allelic pattern associated with this phenotype.

For the same mouse population, Norgard *et al.* (2008) estimated the heritability of long bone length to be ∼0.9, of which we could explain almost half by dimension 1 of the multivariate mapping. The lack of the remaining heritability is due likely to the limited number of SNPs and incomplete linkage disequilibrium between causal variants and genotyped SNPs (Yang *et al.* 2010).

Estimates of explained variance based on all scored loci tend to be too high because of massive overfitting and, hence, may not reflect actual prediction accuracy (Makowsky *et al.* 2011; Gianola *et al.* 2014). After applying leave-one-out cross-validation, dimension 1 (limb length) still showed an explained variance of 0.36 and dimension 3 of 0.17. The reduction of the genetic variables to the first 80 PCs together with the identification of the relevant predictor variable (the genetic latent variable) thus prevented severe overfitting. Our estimates include both additive and dominance effects, but most of the explained variance was due to additive gene effects (fractions of phenotypic variance explained by additive effects were 0.40, 0.05, and 0.23 for the three dimensions).

The variance of body/organ weight (second latent variable) that was explained by the genetic latent variable was relatively small (7% of phenotypic variance), which is somewhat surprising since the two parental strains were selected for small and large body size, respectively. Apart from substantial environmental variance, this may result from the considerable sex interactions identified by Fawcett *et al.* (2008), which we did not include in our analysis. Note also that dimension 2 covers the genetic effects on body/organ weight, *independent* of the effects on limb length, so some of the QTL effects on weight might have been captured by a more general size factor with the limb length. Higher estimates of explained variance of body weight in earlier genome-wide studies likely resulted from overfitting the genotype–phenotype relationship. For example, Kramer *et al.* (1998) explained 47% of phenotypic variance in body weight by a multiple regression on the additive and dominance scores of all scored SNPs. We could reproduce this result with the present sample, but when applying a leave-one-out cross-validation, this fraction dropped to ∼3%. This severe overfitting by the multiple regression is not surprising, given that we found only a single allele combination to be considerably (and presumably causally) associated with body/organ weight. The 705 remaining combinations of genotype scores inflated the “explained variance” by random associations with body weight (note that we had additive and dominance scores for 353 loci and hence 706 independent linear combinations of scores).

Multivariate genotype–phenotype mapping is primarily an exploratory method for investigating the multivariate structure of genotype–phenotype association. However, it can provide crucial information for gene identification. In our mouse data, for instance, we found three independent dimensions of genotype–phenotype association, each of which could be tested as a whole against the null hypothesis of no association. In fact, this is the number of independent statistically meaningful tests that can be made. The performance of all pairwise tests between genetic and phenotypic variables (7766 for our data) is a misuse of significance testing in an entirely exploratory context (McCloskey and Zilik 2009; Bookstein 2014; Mitteroecker 2015) that does not guarantee repeatable results (Morgan *et al.* 2007). Of course, the biological meaning of a hypothesis about the complete lack of genotype–phenotype association remains doubtful nonetheless. Dimensions 1 and 3 were highly statistically significant as a whole in our data; dimension 2 was convincing as a pattern but not as clearly significant as the other dimensions because of the large fraction of environmental variance ( Table 2). Identification of important alleles should be based on the effect sizes (the genetic coefficients), unless more specific prior hypotheses about gene effects existed. The fourth and all subsequent dimensions did not differ significantly from a random association ().

For complex phenotypes measured by multiple variables, multivariate genotype–phenotype mapping should precede any other mapping technique to identify the number of independent dimensions of genotype–phenotype association. In particular, this applies to variables that do not bear biological meaning one by one, such as in modern morphometrics and image analysis, but also to gene expression profiles and similar “big data.” Multivariate genotype–phenotype mapping identifies the phenotypes (linear combinations of measurements) under strong genetic control that are worth considering for further genetic analysis. In addition to the allelic effects estimated by multivariate mapping, other measures of allelic effects or LOD scores can be computed by more classic methods for the identified phenotypic latent variables.

We presented four different variants of multivariate genotype–phenotype mapping, which maximize different measures of association: (i) genetic effect, (ii) genetic variance, (iii) heritability, and (iv) the covariance between genetic and phenotypic latent variables (Table 1). The choice among them depends on the scientific question and the kind of phenotypic variables. While the genetic effect may be of interest in certain medical studies, additive genetic variance is central to many evolutionary studies and breeding experiments. Maximizing heritability tends to be the most unstable approach because it maximizes genetic variance and minimizes environmental variance at the same time. Approach iv, the maximization of covariance, has no correspondence in quantitative genetics but it is computationally simple and avoids severe overfitting without prior variable reduction (Martens and Naes 1989). The genetic coefficients in approaches i–iii represent partial effects (*i.e.*, effects conditional on the other loci), whereas the coefficients in approach iv do not depend on the other loci in this way. Approach iv thus offers an alternative to the other approaches when computational simplicity is preferred over interpretability or when partial coefficients should be avoided.

In addition to purely additive or average gene effects, nonadditive effects can be incorporated into the analysis by adding variables representing dominance or epistasis (pairwise or higher-order interaction terms) to the genetic predictors. Accordingly, the genetic latent variables can be interpreted either as breeding values or as genotypic values. Multivariate genotype–phenotype mapping can be applied to crosses of inbred strains as well as to natural populations. Genetic similarity and common ancestry can be accounted for by generalized least-squares variants (see *Appendix B*, *Generalized least squares*). In addition to genetic variables, covariates such as environmental variables can be included in the predictors as well. In approaches i–iii, the resulting coefficients of the gene effects are then conditional on these covariates. The presented methods make no distributional assumptions and do not require linear relationships between genetic and phenotypic latent variables or covariates. However, only if all relationships are linear (and, hence, the variables jointly normally distributed), uncorrelatedness implies actual independence. The interpretation of the singular values of **F**, **G**, and **H** as genetic effect, genetic variance, and heritability, respectively, is exact only for randomly mating populations in Hardy–Weinberg equilibrium. The more a population deviates from equilibrium, the more the singular values may deviate from these genetic quantities.

For a single phenotypic variable only, approaches i–iii lead to the same genetic latent variables and the same genetic coefficients, which are the regression coefficients of the phenotype on the loci. The corresponding linear combination of loci maximizes all three measures of genotype–phenotype association: genetic effect, genetic variance, and heritability. In the classic genetic literature, this is also known as the “selection index” (Smith 1936; Hazel 1943). The three approaches can thus be construed as three different generalizations of multiple regression to many phenotypic traits; the genetic coefficients (the elements of the vectors ) in all three approaches equal the multiple-regression coefficients of the corresponding *i*th phenotypic latent variable on the loci.

This property allows one to rotate the phenotypic latent variables to increase their biological interpretability, as in exploratory factor analysis, and to estimate the corresponding genetic effects. For the mouse data, the first 3 dimensions of phenotypic latent variables constituted a 3-dimensional subspace of the 11-dimensional phenotype space, which contained almost all of the genetic variation in the data. The first latent variable was overall limb length, and the third one was a contrast between distal and proximal limb bone lengths. Hence, femur and humerus, as well as tibia and ulna, were highly correlated and clustered in Figure 5. Alternative latent variables would thus be proximal limb length () and distal limb length (). The corresponding genetic latent variables can be computed by multiple regression of the new latent variable on the loci.

In most modern genetic data sets *p* clearly exceeds *n*, which challenges least-squares methods such as approaches i–iii. In *Appendix B*, *Computational properties* we show how they can be computed based on generalized inverses or matrix regularizations. Approach iv—the partial least-squares analysis—does not require the inversion of a matrix and can also be applied to collinear genetic variables and when Clearly, there is much room for improvement, such as an implementation in a Bayesian framework and the application of other penalized or BLUE/BLUP methods (*e.g.*, Meuwissen *et al.* 2001; Lopes and West 2004; de Los Campos *et al.* 2013; Zhu *et al.* 2014). The use of information measures, such as the application of PLS to Kullback–Leibler divergences by Bowman (2013), may allow for the application of the presented approaches to a wide range of heterogeneous variables.

In our 11-dimensional phenotype space, only 3 dimensions had considerable genetic variation, but the majority of genetic variation was even concentrated in a single dimension. In this sense, the genotype–phenotype map in this population is of surprisingly “low dimension” (even if the metaphor of dimensionality does not uniquely translate into an integer or a real number). In a preliminary analysis, we found similar results for the F_{9} and F_{10} generations of the same mouse cross. At least in part, this low-dimensional genotype–phenotype map resulted from the intercross of two inbred populations. It remains to be investigated to what degree it is also characteristic of outbred populations. Current studies of phenotypic and genetic variance–covariance patterns provide inconsistent results in this regard (*e.g.*, Kirkpatrick and Lofsvold 1992; Mezey and Houle 2005; Hine and Blows 2006; Pavlicev *et al.* 2009). Hallgrimsson and Lieberman (2008) speculated that a low-dimensional pattern of phenotypic variation is a general phenomenon that results from the “funneling” of the vast amount of genetic variation by a few central developmental pathways and morphogenetic processes. This would massively bias and constrain a population’s phenotypic response to natural or artificial selection and generate a broad heterogeneity of genetic responses within a single selection scenario. If such funneling processes exist, multivariate genotype–phenotype mapping can help to identify these central pathways.

## Acknowledgments

We are very grateful to Fred L. Bookstein, Hans Nemeschkal, and Günter Wagner, as well as to three anonymous reviewers and William Valdar for very helpful comments and discussion. M.P. acknowledges support from the March of Dimes Prematurity Research Center Ohio Collaborative (grant 22-FY14-470).

## Appendices

### Appendix A: Derivation and Proofs

#### Univariate genotype–phenotype association

To show how—for an idealized, randomly mating population—the singular value decompositions of the matrices **F**, **G**, and **H** introduced in (2), (3), and (4), respectively, lead to latent variables with maximal genotype–phenotype association, we first review the classic measures of association in the simple case of one diploid locus with two alleles affecting one quantitative phenotypic trait. For a sample of *n* specimens, let the vector **x** contain the additive genotype scores (1 for heterozygotes and 0 or 2 for the two possible homozygotes) and the vector **y** contain the corresponding phenotypes. Consider further the vector **d** containing the dominance scores (0 for homozygotes and 1 for heterozygotes) and the matrix For notational convenience, let **x**, **y**, and **d** be mean centered so that where the superscript ′ indicates the transpose operation and **1** a vector of 1’s.

In a sample of measured individuals, the additive and dominance effects, *a* and *d*, are numerically identical to the regression coefficients of the multiple regression of the phenotype **y** on both **x** and **d**: (A1)By contrast, the average effect, *α*, of an allele substitution equals the bivariate regression slope of the phenotype **y** on the additive genotype scores **x**: (A2)The sum of the average effects of all alleles constitutes the breeding value of an individual.

The additive genetic variance, of **y** owing to the variance in **x** equals (Falconer and Mackay 1996; Roff 1997). For a population in Hardy–Weinberg equilibrium, and thus the additive genetic variance can be expressed as (A3)where is **x** standardized to unit variance. The dominance variance, is equal to

The additive genetic variance of **y** due to **x**, expressed as a fraction of the total variance of **y**, is the narrow-sense heritability of **y** resulting from variation in the studied locus. It can be expressed as the squared correlation between the locus and the phenotype, (A4)where is the phenotype standardized to unit variance. The broad-sense heritability is the sum of additive and dominance variance as a fraction of phenotypic variance, which equals the squared multiple correlation coefficient resulting from the regression of **y** on **Z**.

All these parameters represent different aspects of the genotype–phenotype relationship. While *a* and *d* are properties of the genotype alone, *α* and are population properties that represent the potential to respond to natural or artificial selection. In contrast to *a* and *d*, the average effect *α* depends on the allele frequencies and in a population: The heritability depends on both genetic and nongenetic variation in the population (see also Hansen *et al.* 2011).

#### Multivariate genotype–phenotype association

For multiple loci and multiple phenotypic variables, we seek the combined effect of the alleles on a phenotype composed of the measured variables In a cross of two inbred lines, each locus is represented by one variable for the additive genotype scores and one for the dominance scores. In natural populations, where each locus can have multiple alleles, each allele at each locus is represented by a separate variable containing the number of this allele (0, 1, 2) at the locus and one variable for the dominance scores. In the following notation, the variables **x**_{i}, **y**_{i} are again assumed to be mean centered. Let the effects of the alleles on the phenotype be denoted by the vector and the weightings of the measured variables that determine the phenotype by the vector Then the two latent variables are given by the linear combinations and (Figure 1). The coefficient vectors are chosen to maximize the association between the corresponding latent variables,where *A* represents one of the above association functions (genetic effect, genetic variance, heritability) between the genetic and phenotypic latent variables. In addition to this pair of latent variables, there might be further pairs of latent variables with effects and independent of the previous ones, that together account for the observed genotype–phenotype association.

#### Maximizing genotype–phenotype association via SVD

The association functions (A1)–(A4) extend naturally from a single locus and a single trait to a linear combination of loci and a linear combination of phenotypic variables.

##### Genetic effect (i):

Under the constraint the regression slope of on conditional on all other linear combinations of **X**, is maximized by the first pair of singular vectors of the matrix of multiple-regression coefficients If **X** contains the *p* additive scores only, the singular value represents the average effect associated with the allele combination and the phenotype specified by the singular vectors. Whereas if **X** comprises both additive and dominance scores ( in total), the singular value is the sum of additive and dominance effects.

To prove this, consider the matrix of regression coefficients for the linear combinations and where **A** and **B** are orthonormal matrices containing the vectors and respectively: (A5)(Note that this equation holds only under the constraint that the vectors are mutually orthogonal so that the matrix **A** is orthonormal and ) The right part of Equation A5 is a classic singular value problem (*e.g.*, Mardia *et al.* 1979). If **A** is equal to the matrix **U** of left singular vectors of **F** and **B** is equal to the matrix **V** of right singular vectors, then is the diagonal matrix of singular values of **F**. These singular values are equal the regression slopes of the phenotypic latent variables on the corresponding genetic latent variables. The pair of singular vectors associated with the largest singular value determines the pair of genetic and phenotypic latent variables with maximal regression slope (*i.e.*, maximal genetic effect),

##### Genetic variance (ii):

The variance in the phenotype explained by the allele combination is given by (compare Equation A3). It is maximized by the vectors and where and are the first left and right singular vectors of The matrix **G** is equal to the covariance matrix between **Y** and **X** after **X** is transformed to a spherical distribution. If **X** contains the additive scores only, is the additive genetic variance associated with the corresponding allele combination and phenotype, whereas if **X** contains both additive and dominance scores, it represents the total genetic variance.

To show this, express the maximization of as the maximization of subject to Write then the maximization is of subject to This is well known to be solved by the first pair of singular vectors of where and with the maximal variance equal to (*cf*. Mardia *et al.* 1979; p. 284).

##### Heritability (iii):

The squared correlation coefficient is maximized by the vectors and where and are the left and right singular vectors of the matrix the covariance matrix between the spherized **Y** and spherized **X**. Depending on whether **X** contains the additive scores or both additive and dominance scores, the squared singular value (squared canonical correlation) is the narrow-sense or broad-sense heritability, respectively. The proof follows by extending the one of approach ii; it equals the standard proof of canonical correlation analysis (*e.g.*, in Mardia *et al.* 1979).

For each of the three approaches, further dimensions (pairs of latent variables) can be extracted by the subsequent pairs of singular vectors of the corresponding association matrix:(i)(ii)(iii)In approach i, the genetic effects (the vectors ) are orthogonal because they are the singular vectors of a real matrix. In the maximizations of (ii) genetic variance and (iii) heritability, by contrast, the genetic latent variables are uncorrelated. This can be shown by expressing the covariance matrix of the genetic latent variables as where the orthonormal matrix **U** contains the singular vectors and **I** is the identity matrix. The phenotypic effects are orthogonal in approaches i and ii, whereas the phenotypic latent variables are uncorrelated in approach iii.

In all three approaches, the corresponding association is maximized conditional on all other linear combinations of alleles, including the other dimensions This implies that where is the regression slope of the *i*th phenotypic latent variable on the *i*th genetic latent variable and is the regression slope conditional on another linear combination for any vector The matrix of regression coefficients of the phenotypic latent variables on the genetic latent variables thus is diagonal.

The proof of these properties for approach i follows from the fact that the matrix of regression coefficients for the latent variables, and can be written as (see Equation A5), where **Λ** is the diagonal matrix of singular values of **F**. For approach ii, and where **U** and **V** are the orthonormal matrices of left and right singular vectors of **G**. The matrix of regression coefficients for the latent variables can thus be written as (A6)where the diagonal matrix **Λ** now contains the singular values of **G**. The proof for approach iii can be constructed similarly. Furthermore, in approaches ii and iii the genetic latent variables are correlated only with the corresponding phenotypic latent variable but not with any of the other latent variables: for This can be shown by expressing the cross-covariance matrix of the latent variables as in approach ii and as in approach iii, which both are diagonal.

##### Partial least-squares analysis

The fourth approach, PLS, maximizes the covariance between the two linear combinations The unit vectors maximizing this covariance can be computed as the first pair of singular vectors of the cross-covariance matrix The proof of this classic singular value problem is given, *e.g.*, in Mardia *et al.* (1979); see also Sampson *et al.* (1989). Subsequent pairs of singular vectors yield further genetic and phenotypic dimensions that are mutually orthogonal. Furthermore, for and for where is the *i*th singular value of

The scaled genetic coefficients (the elements of the scaled singular vectors are equal to the covariances between the corresponding locus and the phenotypic latent variable, without conditioning on the other loci as in approaches i–iii. When the genetic variables are standardized to unit variance through division by their standard deviation, the squared scaled genetic coefficients equal the explained variance of the phenotypic latent variable owing to the corresponding locus considered separately, *i.e.*, without conditioning on the other loci (compare Equation A3).

### Appendix B: Properties

#### Geometric properties

The maximizations of genetic effect and of covariance in approaches i and iv require a constraint on the length of **a** and **b**. Because they are the singular vectors they are computed to have a 2-norm of 1. When the variables can be equipped with a meaningful Euclidean metric, this constraint translates into an interpretable notion of total genetic and phenotypic effects. For the genetic variables, this choice of constraint is not particularly obvious as it implies that an allele with an effect of 1 is equivalent in magnitude to two alleles with effects of each. It may seem more intuitive that two alleles with effects of are equivalent to a single allele with an effect of 1. This latter choice would impose a constraint on the sum of the absolute values of the elements of **a** (the 1-norm), not on the sum of the squared elements. Regression approaches with constraints other than the 2-norm have been proposed [*e.g.*, the Lasso technique (Tibshirani 1996)], but this generalization goes beyond the scope of the present article.

The invariance to the length of **a** in approach ii can also be seen from Equation A3, which expresses the additive genetic variance of a single trait as its squared covariance with the additive genotype scores scaled to unit variance. The scale of **x** is removed through dividing by its standard deviation. In the multivariate context, the maximal genetic variance, can be found by maximizing where is the matrix of genetic variables transformed so that every linear combination has unit variance: for any unit vector **c**. This transformation is achieved by multiplying **X** by the inverse square root of its covariance matrix: The vectors **a** and **b** can thus be found by the singular vectors of after has been transformed back into the original coordinate system. The singular values of this matrix (the genetic variances) as well as the right singular vectors determining the phenotypes are invariant to affine transformations of **X**. This can be shown when considering that the right singular vectors and squared singular values of **G** are the eigenvectors and eigenvalues of (*cf*. Equation B1). Let be a linear transformation of **X**, where **T** is a full-rank matrix, and It follows from this property that the maximal genetic variances (singular values of **G**) as well as the phenotypes that show these maximal genetic variances (right singular vectors) remain unchanged by linear transformations of the genetic variables; thus they also do not depend on the variances and covariances of the genetic variables, that is, on genetic variance and linkage disequilibrium. The same property holds for approach iii. Approaches i and iv are not invariant to transformations of **X**, implying that the genetic variances and covariances need to be interpretable for computing the maximal genetic effects. Only approach iii is invariant to affine transformation of the phenotypic variables **Y**; all other approaches require meaningful phenotypic variances and covariances as well as commensurate units.

#### Computational properties

For typical genetic data, approaches i–iii are not computable by the presented least-squares methods because of collinearities between loci or because The covariance matrix is singular and its inverse or inverse square root cannot be computed without the use of a pseudoinverse or of regularization techniques. The Moore–Penrose pseudoinverse of a matrix **M** is where **Q** is the matrix of eigenvectors of **M** and is a diagonal matrix with the reciprocal of the *m* largest nonzero eigenvalues in the diagonal. The remaining eigenvalues are set to 0. This is equivalent to reducing the data to the first *m* principal components for the inversion, discarding all subsequent principal components with small or zero variance. The partial coefficients resulting from such an approach are not conditional on all other variables, but only on the major patterns of multivariate variation (discarding rare alleles). In the simplest form of Tikhonov regularization, also referred to as ridge regression, is replaced by where **I** is the identity matrix and *γ* is a positive real. The larger *γ* is, the more are components with low variance downweighted in the matrix inverse; *i.e.*, rare alleles or less variable allele combinations are downweighted relative to more variable alleles or allele combinations.

The results of approaches i–iii depend on the number of selected components or on *γ* and require a careful decision. Typically, after adding the first few components that cover the relevant signals, adding further components has little effect, until the number of components becomes too large and the increasing noise leads to unstable results. Adding further components may still increase the explained phenotypic variance, but this is due to overfitting; the genetic coefficients may not be interpretable. Exploring different numbers of components underlying the pseudoinverse (or different values of *γ*) thus often leads to a range of stable components (or a stable range of *γ*) that lead to similar and equally interpretable results. A cross-validation approach can help to find the optimal number of principal components or the optimal *γ* for the data set. Approach iv—the partial least-squares analysis—involves no matrix inverse and can also be computed for collinear loci and if overfitting is less a problem in this approach (Martens and Naes 1989). It is thus useful to compare the results of approaches i–iii to that of approach iv. Bayesian approaches and numerous other penalized methods offer promising alternatives to the presented least-squares methods (Meuwissen *et al.* 2001; Lopes and West 2004; de Los Campos *et al.* 2013; Zhu *et al.* 2014).

If the number of genetic variables *p* or the number of phenotypic variables *q* is very large, the singular value decomposition of the association matrix can be computationally demanding. Here one can make use of the property that the left singular vectors of a matrix **M** are equal to the eigenvectors of and the right singular vectors are the eigenvectors of Thus, if or one can compute either or as the eigenvectors of the smaller matrix product. Since and the other singular vectors can be obtained by premultiplication with **M**. In approach ii, for instance, the vectors are given by the eigenvectors of (B1)Note that the eigenvalues of (B1) are equal to the squared singular values of **G**.

If both *p* and *q* are very large and if only the first few dimensions need to be computed, the singular vectors can be computed more effectively via an iterative approach. Start with any vector and estimate as scaled to unit vector length. In the next step, is estimated as again scaled to unit vector length. These steps are repeated until convergence, which is usually reached fast. The singular value equals To compute the next pair of singular vectors let and repeat the iterative approach with instead of and similarly for subsequent dimensions.

#### Generalized least squares

In a sample with a family structure, the presented least-squares estimates are unbiased but the standard errors of the parameters may be inflated. This can be addressed by generalized least squares. Let the matrix **Ω** contain measures of expected or realized genetic relatedness between pairs of individuals (*e.g.*, Hayes *et al.* 2009). Then the maximization in the four approaches is based on the following matrices:

#### The genetic variance–covariance structure

Because in approach ii—the maximization of genetic variance—the vectors constitute an orthonormal basis of the phenotype space with associated genetic variances the genetic variance–covariance matrix of the measured phenotypic variables can be computed as where **B** contains the phenotypic coefficient vectors (the left singular vectors of **G**), and is a diagonal matrix of the squared singular values of Because of (B1), the genetic variance–covariance matrix can also be computed directly as (B2)This can also be seen when considering the phenotypic predictions from the genetic variables, of which is the variance–covariance matrix. Note that when **X** contains the additive genotype scores only, (B2) is the additive genetic variance–covariance matrix, whereas if **X** contains additive and dominance scores, (B2) is the total genetic covariance matrix.

### Appendix C: Alternative Analyses

Here we present the results of the other three approaches applied to the same mouse data as in the main text. All three approaches lead to similar scree plots as in approach ii: The first dimension clearly dominates the genotype–phenotype relationship (Figure C1). In all approaches, the first dimension represents limb length (Figure C2) and, hence, also the genetic coefficients are highly consistent (not shown).

Approach iii, the maximization of heritability, is most unstable and requires a dimension reduction of the phenotypic variables. In this analysis we used the first four principal components of the 11 measurements. By contrast, the partial least-squares analysis in approach iv does not require variable reduction or regularization for the genetic and phenotypic variables.

The second dimension largely reflects body and organ weight in all four approaches (Figure C2), and also the genetic coefficients are very similar (not shown). The third dimension differs among the approaches. Only in approach iv the third dimension reflects a contrast between distal and proximal long bones as in approach ii.

## Footnotes

*Communicating editor: W. Valdar*Supplemental material is available online at www.genetics.org/lookup/suppl/doi:10.1534/genetics.115.181339/-/DC1.

- Received July 29, 2015.
- Accepted February 15, 2016.

- Copyright © 2016 by the Genetics Society of America