I'm working on a biological question, with species data derived from an external database, which has multiple response and predictor variables. As a result, I want to do multivariate regression across a phylogeny to empirically test if my response variables are significantly different in respect to my predictors.
Please refer to source  and it's citations for your own investigation of this process.
I know how to do multiple regression via pGLS, but the R package  only mentions predictors and response. Furthermore, another source  discusses how multivariate regression though pGLS in R can be done, but requires one to transform the data under a Brownian motion model. (Edit: It seems that  is a solution...so I'm looking the process).
The vignette for the pGLS package (pdf)
D.C. Adams. 2014. A Method for Assessing Phylogenetic Least Squares Models for Shape and Other High-Dimensional Multivariate Data. Evolution. 68:9 2675-2688. doi: 10.1111/evo.12463
Revell, L. J. (2010), Phylogenetic signal and linear regression on species data. Methods in Ecology and Evolution, 1: 319–329. doi: 10.1111/j.2041-210X.2010.00044.x