We consider a $d$-dimensional branching particle system in a random environment. Suppose that the initial measures converge weakly to a measure with bounded density. Under the Mytnik-Sturm branching mechanism, we prove that the corresponding empirical measure $X_t^n$ converges weakly in the Skorohod space $D([0,T];M_F(\mathbb{R}^d))$ and the limit has a density $u_t(x)$, where $M_F(\mathbb{R}^d)$ is the space of finite measures on $\mathbb{R}^d$. We also derive a stochastic partial differential equation $u_t(x)$ satisfies. By using the techniques of Malliavin calculus, we prove that $u_t(x)$ is jointly Hölder continuous in time with exponent $\frac{1}{2}-ε$ and in space with exponent $1-ε$ for any $ε>0$.