Given \(p\ge 2\) and a map \(g : B^n(0,1)\rightarrow S_n^{++}\) , where \(S_n^{++}\) is the group of positively definite matrices, we study critical points of the following functional: \(v\in W^{1,p}\left( B^n(0,1);\mathbb {R}^N \right) \mapsto \int _{B^n(0,1)} |\nabla v|^p_g\, d\textrm{vol}_g = \int _{B^n(0,1)} \left( g^{\alpha \beta }(x) \left\langle \partial _\alpha v(x),\partial _\beta v(x) \right\rangle \right) ^{\frac{p}{2}}\, \sqrt{\det g(x)}\, dx.\) We show that if g is uniformly close to a constant matrix, then v is locally Hölder-continuous. If g is Hölder-continuous, we show that \(\nabla v\) is locally Hölder-continuous. As an application, we prove that any Hölder-continuous solution to \(|\Delta _{g,p}u|\lesssim |\nabla u|^p_g\) satisfies additional regularity properties depending on the regularity of g. In the case \(p=n\) , only the continuity is assumed a priori.