In biomedical research, multivariate functional data are frequently encountered. Majority of the existing approaches for functional data analysis focus on univariate functional data and the methodology for multivariate functional data is far less studied. Particularly, the problem of investigating covariate effects on multivariate functional data has received little attention. In this research, we propose a fully Bayesian latent factor regression for studying covariate effects on multivariate functional data. The proposed model obtains a low-dimensional representation of multivariate functional data through basis expansions for splines and factor analysis for the basis coefficients. Then, the latent factors specific to each functional outcome are regressed onto covariates accounting for residual correlations among multiple outcomes. The assessment of covariate effects is conducted based on the marginal inclusion probability for each covariate, which is calculated a posteriori by assigning a stochastic search variable selection (SSVS) prior to the regression coefficients. To better control for the false discovery rate, we propose a multivariate SSVS prior that allows for a set of coefficients to be zero simultaneously. We illustrate the proposed method through a simulation study and an application to the air pollution data collected for 13 cities in China.