PHENIX: An R package to estimate a size-controlled phenotypic integration index.

PREMISE OF THE STUDY: Organisms usually show intercorrelations between all or some of their components leading to phenotypic integration, which may have deep consequences on the evolution of phenotypes. One of the main difficulties with phenotypic integration studies is how to correct the integration measures for size. This has been considered a challenging task. In this paper, we introduce an R package (PHENIX: PHENotypic Integration indeX), in which we provide functions to estimate a size-controlled phenotypic integration index, a bootstrapping method to calculate confidence intervals, and a randomization method to simulate null distributions and test the statistical significance of the integration. METHODS AND
RESULTS: PHENIX is an open source package written in R. As usual for R packages, the manual and sample data are available at: http://cran.r-project.org/web/packages/PHENIX/index.html. Functions included in this package easily estimate phenotypic integration by controlling a third variable (e.g., the size of the studied organ).
CONCLUSIONS: PHENIX helps to estimate and test the statistical significance of the magnitude of integration using one of the most-used methodological approaches, while taking size into account.

Organisms consist of complex phenotypes that usually show intercorrelations between all or some of their components leading to phenotypic integration and modularity (Berg, 1960; Pigliucci, 2003). Interconnections between phenotypic traits reflect webs of developmental, physiological, and genetic interactions behind the genotype-to-phenotype map (Murren, 2002; Armbruster et al., 2014). If a suite of traits is integrated, important questions arise. For instance, is the observed integration the result of natural selection or does it reflect developmental or genetic relationships among traits limiting and channeling its evolutionary course (Pigliucci, 2003; Hansen and Houle, 2004)? The integration of the phenotypes then has functional implications that may shape the evolution of phenotypes.

In plants, most of the research on phenotypic integration has been focused on flowers. Flowers are complex structures composed of different organs that together work toward the same aim: offspring production. These floral organs may show different degrees of phenotypic integration (Ordano et al., 2008; Harder, 2009). Several factors seem to cause floral integration. First, floral phenotypic integration may result from the effects of genetic and/or developmental constraints (Ashman and Majetic, 2006; Bissell and Diggle, 2010). Second, stabilizing selection imposed by pollinators on different floral organs may drive a decrease of phenotypic variation and an increase of phenotypic integration (Armbruster, 1991; Pérez-Barrales et al., 2007; Rosas-Guerrero et al., 2011). Whatever the cause, floral integration may have deep consequences for floral trait evolution, such as the floral traits involved in the accurate placement of pollen grains on pollinators or stigmas (Armbruster, 1991; Pérez-Barrales et al., 2007; Pérez et al., 2007; Bissell and Diggle, 2010; Ferrero et al., 2011; Rosas-Guerrero et al., 2011).

Although phenotypic integration is an important concept to understand the evolutionary process (e.g., Pigliucci and Preston, 2004), it is also a complex concept (Magwene, 2001; Pavlicev et al., 2009; Haber, 2011; Armbruster et al., 2014) and difficult to distill into a single index. For instance, integration values might be, in part, a consequence of resource acquisition variability between individuals. This effect is more evident when the traits analyzed are resources allocated to different components of a given organ or individual. When resource availability is not controlled, larger individuals may have larger components, driving a correlation between components due, in part, to variation in resource availability. This correlation may obscure others that are unrelated to size, such as genetic correlations. To disentangle whether the levels of integration observed in a population are due to genetic constraints or to resource availability, Torices and Méndez (2014) proposed a size-controlled integration index using a modification of the widely used Wagner (1984) index, in which the size of the studied structure, as a proxy of resource availability, is taken into account.

In this paper, we introduce the PHENIX package (PHENotypic Integration indeX) available for the open sourced statistical software R (R Core Team, 2014). In this package, we provide functions to estimate the phenotypic integration index proposed by Wagner (1984), where size cannot be considered, and the new size-controlled index (Torices and Méndez, 2014). In addition, a bootstrapping method to calculate confidence intervals and a randomization method to test the statistical significance of the integration are available for both indices. Finally, we provide an example using data on resource allocation to floral components from Méndez and Traveset (2003) to show its applicability.

METHODS AND RESULTS

PHENIX estimates the phenotypic integration index proposed by Wagner (1984) and includes a modification to this index proposed by Torices and Méndez (2014) where the effect of a third variable under study (e.g., size) can be taken into account. In the phenotypic integration index proposed by Wagner (1984) (PINT hereafter), the magnitude of phenotypic integration among a set of traits is quantified by the variance of the eigenvalues (λ) of the correlation matrix between traits (Wagner, 1984; Cheverud et al., 1989). For the size-controlled index (PINTsc hereafter), the correlation matrix is replaced by a partial-correlation matrix (Torices and Méndez, 2014). With the partial-correlation approach, the correlation structure between all traits is assessed after controlling by a third variable on the correlation between the studied traits.

PINT value can only be directly comparable between data sets containing the same number of traits (N) and individuals (n) because the expected random covariation among traits depends on these values (specifically, it is determined by (N − 1) / n; Wagner, 1984; Cheverud et al., 1989; Pavlicev et al., 2009). Thus, it is common to correct PINT by subtracting the expected amount of PINT produced by random covariation. In the PHENIX package, the function ‘pintsc’ estimates (i) the observed PINTsc (Eq. 1), (ii) the corrected PINTsc in which the expected amount of integration produced by random covariation is removed (Eq. 2), and (iii) a relative value in which the magnitude of the phenotypic integration is also expressed as a percentage of the maximum possible value, by scaling PINTsc values according to the maximal eigenvalue variance determined by the number of measured traits minus 1 (Pavlicev et al., 2009; Eq. 3). The PINTsc confidence intervals are calculated by bootstrapping using the function ‘pintsc.boot’. The amount of random covariation between traits in partial-correlation matrices (Eq. 2) was assumed to be the same as in correlation matrices. However, the effect of the structure of the partial-correlation matrices on the variance of the eigenvalues requires further scrutiny to confirm that this correction can be applied to PINTsc. In addition, the functions ‘pint’ and ‘pint.boot’ also estimate the PINT index, without control of a third variable, and its confidence intervals, respectively. Eqs. 1, 2, and 3 are the same for both PINT and PINTsc with the only difference being that λ represents the variance of the eigenvalues for a correlation matrix in the former and the variance of the eigenvalues of a partial-correlation matrix in the latter.

The statistical significance of both integration indices, PINT and PINTsc, can be assessed by means of the functions ‘pint.p’ and ‘pintsc.p’, respectively. These functions generate a random-correlation matrix with a number of rows and columns equal to the number of traits in the input data set. Diagonal elements in this square symmetric matrix are set to 1, whereas off-diagonal values are randomly selected from a Pearson product moment correlation coefficient distribution (as implemented in the SuppDists package [Wheeler, 2013]). This distribution is defined between −1 and 1, and its shape may be controlled by the user with the ‘N.Pearson’ parameter. Thus, the probability of drawing a value to fill the simulated correlation matrix will be higher for extreme values (that is, those closer to −1 and 1) if the ‘N.Pearson’ parameter is set to 3 or for central values (that is, those closer to 0) if ‘N.Pearson’ is higher than 4. If ‘N.Pearson=4’, values will be drawn from a uniform distribution (that is, all values between −1 and 1 show the same probability of being sampled). However, Harder (2009) demonstrates that using a uniform distribution overestimates the expected integration index. We selected 15 as a default value for the ‘N.Pearson’ parameter because it generates distributions in line with those observed by Harder (2009).

The resulting simulated random matrix is used to estimate an integration index as implemented in ‘pint’ or ‘pintsc’. The procedure is repeated according to a value defined by the user (1000 times by default) to obtain a set of simulated integration values to be used as a null distribution under the hypothesis of random correlations between every pair of traits. Given a set of traits, the null distribution for PINT and PINTsc is expected to be the same under random correlation among all traits (Appendix S1). Both distributions associated with this method (that is, the one defined by the ‘N.Pearson’ parameter and the resulting null distribution for the integration index) can be visualized using the ‘plot’ option (set to ‘P’ and ‘R’, respectively).

We provide an example of the potential applications of PHENIX using a data set from a study on resource allocation to different flower components in the mostly single-flowered Paeonia cambessedesii (Willk.) Willk. (Paeoniaceae), performed by Méndez and Traveset (2003). This species is a self-compatible, herbaceous perennial plant endemic to the Balearic Islands (Mediterranean Sea). Among other things, Méndez and Traveset (2003) studied the patterns of correlation between flower components using three allocation currencies (dry mass, nitrogen [N], and phosphorus [P]). Dry mass and P allocation to stamens and gynoecium were positively correlated, thus suggesting that resource allocation to floral components could be significantly integrated. Two plausible explanations arise from this pattern: (i) resource allocation to floral components might be correlated because of genetic correlations or (ii) the observed correlation between resource allocation to floral components might also be the result of differences in size between flowers included in the study. That is, larger flowers might also display both larger stamens and larger gynoecia compared to smaller flowers, driving this positive correlation between allocations to different components. To assess this effect, the size of the flower can be used to control for the resource acquisition variation. The comparison between PINT and PINTsc indices (without and with size-controlling, respectively) may help to assess the role of size in the observed magnitude of integration.

We estimated PINT and PINTsc for the P. cambessedesii data set using PHENIX functions and the example data set ‘paeonia’, included in the package:

# to load Mendez and Traveset's (2003) data on dry  mass allocation to floral components:

data(paeonia)

# This data set contains nine columns. The first column  indicates the plant ID and the last one represents the total size of each flower (measured as total dry mass). Thus, floral component variables are columns 2 to 8:

fl.components<-paeonia[,2:8]

# to estimate Wagner's integration index:

PIN<-pint(paeonia[,2:8])

PIN

# to estimate confidence intervals by bootstrapping:

PIN.BOOT<-pint.boot(paeonia[,2:8], replicates=5000)

PIN.BOOT

# to test the statistical significance of the integration:

PIN.P<-pint.p(paeonia[,2:8],n.replicates=5000)

PIN.P$Summary

# to save plot as a pdf file (Fig. 1A):

Fig. 1.

Significance level for the phenotypic integration indices of the dry mass allocation to floral components in Paeonia cambessedesii estimated without (A) and with (B) size-controlled correction. Histograms represent the distribution of each index assuming random correlation between traits. The red dashed line represents the phenotypic integration value observed for the real data set. PINT = phenotypic integration index proposed by Wagner (1984); PINTsc = size-controlled index proposed by Torices and Méndez (2014); RelPINT and RelPINTsc = relative values in which the magnitude of the phenotypic integration is expressed as a percentage of the maximum possible integration value; PINT.c and PINTsc.c = PINT and PINTsc indices corrected by subtracting the expected amount of integration produced by random covariation.

dev.copy2pdf(file="paeonia_pint.p.pdf")

# to estimate the size-controlled integration index  it is mandatory to define the “control” variable. In this data set, the control variable is the total size of the flower:

PINSC<-pintsc(traits=paeonia[,2:8], control=paeonia[,9])

PINSC

# to estimate confidence intervals by bootstrapping:

PINSC.BOOT<-pintsc.boot(traits=paeonia[,2:8],  control=paeonia[,9], replicates=5000)

PINSC.BOOT

# to test the statistical significance of the  integration:

PINSC.P<-pintsc.p(paeonia[,2:8], control=paeonia[,9],  n.replicates=5000)

PINSC.P$Summary

# to save plot as a pdf file (Fig. 1B):

dev.copy2pdf(file="paeonia_pintsc.p.pdf")

# to make a custom output (first row in Table 1)

Table 1.

Degree of phenotypic integration among absolute allocation of resources to floral components for three currencies (dry mass, N, and P) in the single-flowered species Paeonia cambessedesii.

	Non–size-controlled index		Size-controlled index
	Corrected PINT ± SE (99% CI)	Rel PINT	Corrected PINTsc ± SE (99% CI)	Rel PINTsc
Dry mass	0.852 ± 0.096*** (0.228−2.363)	18.96***	0.206 ± 0.047ns (0.193−1.324)	8.19ns
Nitrogen	0.828 ± 0.109** (0.200−2.585)	25.46**	0.477 ± 0.062* (0.288−1.733)	16.68*
Phosphorus	0.625 ± 0.100** (0.075−2.370)	20.40**	0.548 ± 0.060* (0.331−1.831)	18.46*

Note: Corrected PINT = phenotypic integration index proposed by Wagner (1984) corrected by subtracting the expected amount of integration produced by random covariation; Corrected PINTsc = the size-controlled index proposed by Torices and Méndez (2014) corrected by subtracting the expected amount of integration produced by random covariation; Rel PINT and Rel PINTsc = relative values in which the magnitude of the phenotypic integration is expressed as a percentage of the maximum possible integration value; SE = standard error; 99% CI = confidence interval at 99%.

The statistical significance of the integration indices was assessed by randomization: ns = nonsignificant; ***P < 0.001, **P < 0.01, *P < 0.05.

row1<-matrix(c(

paste(PIN$PINT.c," ± ", PIN.BOOT[4,2]," (", PIN. BOOT[5,2], "-", PIN.BOOT[6,2],")",sep=""),

PIN.P$Summary[3,3],

PIN$RelPINT,

PIN.P$Summary[2,3],

paste(PINSC$PINTsc.c," ± ", PINSC.BOOT[4,2],"  (", PINSC.BOOT[5,2], "-", PINSC.BOOT[6,2],")",sep=""),

PINSC.P$Summary[3,3],

PINSC$RelPINT,

PINSC.P$Summary[2,3]

),nrow=1)

colnames(row1)<-c("Corrected PINT ± SE (99% CI)",  "P-value","Rel PINT","P-value","Corrected PINTsc ± SE (99% CI)","P-value","Rel PINTsc","P-value")

row.names(row1)<-"Dry mass"

as.data.frame(row1)

We observed that resources allocated to flower components were phenotypically integrated (Table 1). The magnitude of the integration was significantly higher than the magnitude expected for all currencies (Table 1; Fig. 1). The inclusion of size in the estimation of the magnitude of the integration produced a reduction in the observed integration in all studied currencies (Table 1). Although the confidence intervals of both indices (PINT and PINTsc) overlapped, suggesting that the difference between both indices was not statistically significant (Table 1), the statistical tests indicate that at least for dry mass, the observed integration may be significantly due to differences in size.

The reduction in the value of integration when size was taken into account was more pronounced for dry mass allocation, indicating that a higher part of the observed integration could be produced by differences in size, compared to the observed integration for P allocation, which almost did not vary after having controlled by size (Table 1). In addition, the size-controlled measure of integration for dry mass allocation was not significantly different from a null distribution assuming that there is a random correlation between traits (Fig. 1). Therefore, when size was taken into account, the magnitude of the integration in terms of dry mass allocation was not distinguishable from a random pattern of correlation. This result suggests that, to a certain extent, the integration observed in the resources allocated to floral components in P. cambessedesii could be produced by differences in size. Nevertheless, a significant amount of integration was still observed in the allocation patterns to floral components in terms of N and P, even when the potential effect of size was controlled, thus supporting the assertion that genetic correlations may be leading to positive correlations between floral components, such as allocation to male and female organs.

CONCLUSIONS

PHENIX provides a simple solution to assess the effect of size on the magnitude of phenotypic integration using a modification of one of the most-used indices (the Wagner [1984] index, PINT) to quantify the magnitude of integration. One of the main problems with integration studies is how to correct the phenotypic integration measures for size, which has been considered a “thorny task” (Armbruster et al., 2014). The size-controlled index included in PHENIX (PINTsc) takes into account the effect of size on the correlation matrix. Thus, the comparison between both indices (PINT and PINTsc) will inform us of the effect of size on integration. A priori, we could expect a reduction in the magnitude of integration after controlling by size, although the opposite has also been observed (Torices and Méndez, 2014).

PHENIX does not solve problems associated with how to measure size, but it provides a framework to incorporate size once it has been measured. Which metric measure(s) better describes the size of plant organs can be case specific. However, we propose to measure size in terms of dry weight when possible (e.g., Pérez-Harguindeguy et al., 2013). We recognize that it is not always possible to perform destructive harvests to estimate dry weights. However, we recommend looking for proxies in a subset of individuals, exploring what metric trait (e.g., organ length) or combination of traits (e.g., organ length and organ width) can have a strong linear relationship with the organ’s dry weight to be considered subsequently in the estimation of the integration magnitude.

Although this method was developed to control by size, our functions allow the user to control for variables other than size. Thus, this feature broadens the utility of this method beyond size-correction comparisons. Overall, the package and procedures described in this paper may improve the study of phenotypic integration by providing researchers with a framework to estimate and test the statistical significance of the magnitude of the integration with free statistical software.

Click here for additional data file.