Linear spectral transformations of orthogonal polynomials in the real line, and in particular Geronimus transformations, are extended to orthogonal polynomials depending on several real variables. Multivariate Christoffel-Geronimus-Uvarov formulae for the perturbed orthogonal polynomials and their quasi-tau matrices are found for each perturbation of the original linear functional. These expressions are given in terms of quasi-determinants of bordered truncated block matrices and the 1D Christoffel-Geronimus-Uvarov formulae in terms of quotient of determinants of combinations of the original orthogonal polynomials and their Cauchy transforms, are recovered. A new multispectral Toda hierarchy of nonlinear partial differential equations, for which the multivariate orthogonal polynomials are reductions, is proposed. This new integrable hierachy is associated with non-standard multivariate biorthogonality. Wave and Baker functions, linear equations, Lax and Zakharov-Shabat equations, KP type equations, appropriate reductions, Darboux/linear spectral transformations, and bilinear equations involving linear spectral transformations are presented. Finally, the paper includes an Appendix devoted to multivariate Uvarov transformations. Particular attention is paid to 0D-Uvarov perturbations and also to the 1D-Uvarov perturbations, which require of the theory of Fredholm integral equations.