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 (apps.1400104_s1.pdf)). 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:

## 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.*^{a}

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.

## LITERATURE CITED

*Dalechampia scandens*.

*Evolution*45: 1229–1244. Google Scholar

*Philosophical Transactions of the Royal Society of London*.

*Series B, Biological Sciences*369: 20130245. Google Scholar

*Heredity*96: 343–352. Google Scholar

*Evolution*14: 171–180. Google Scholar

*Nicotiana*: Quantitative genetic and comparative phenotypic approaches to floral integration.

*Journal of Evolutionary Biology*23: 1744–1758. Google Scholar

*Systematic Zoology*38: 201–213. Google Scholar

*Lithodora*and

*Glandora*(Boraginaceae).

*Plant Biology*13: 7–18. Google Scholar

*Evolutionary Biology*38: 476–488. Google Scholar

*In*M. Pigliucci and K. Preston [eds.], Phenotypic integration: Studying the ecology and evolution of complex phenotypes, 130–154. Oxford University Press, New York, New York, USA. Google Scholar

*New Phytologist*183: 247–248. Google Scholar

*Evolution*55: 1734–1745. Google Scholar

*Oecologia*137: 69–75. Google Scholar

*Plant Species Biology*17: 89–99. Google Scholar

*New Phytologist*179: 1183–1192. Google Scholar

*Evolutionary Biology*36: 157–170. Google Scholar

*Schizanthus*(Solanaceae): Does pollination truly integrate corolla traits?

*Journal of Evolutionary Biology*20: 1730–1738. Google Scholar

*Narcissus papyraceus*(Amaryllidaceae).

*Oikos*116: 1904–1918. Google Scholar

*Australian Journal of Botany*61: 167–234. Google Scholar

*Ecology Letters*6: 265–272. Google Scholar

*Ipomoea*.

*Evolution*65: 350–364. Google Scholar

*International Journal of Plant Sciences*175: 713–723. Google Scholar

*Journal of Mathematical Biology*21: 77–95. Google Scholar

## Notes

[1] This work was partially supported by a postdoctoral fellowship (BVA 2010–0375) from the Spanish Ministry of Education awarded to R.T. and by the European Regional Development Fund under the NSR Framework awarded to A.J.M.-P. (Program ON.2—O Novo Norte 2007/2013). The authors especially thank Marcos Méndez for kindly providing data on floral allocation for testing the functions included in PHENIX and for comments on manuscript drafts and Phillippa Bennett for correcting the English.