It is non-local because the time evolution of the wave function depends instantaneously on parts of the field which are arbitrary far away. In the case of a linearization using the square root you get.
\begin{equation}
i \hbar \frac{\partial}{\partial t}\, \psi ~~=~~ \tilde{H}\,\psi ~~=\ \sqrt{\mathbf{\tilde{p}}^2 c^2 + m^2c^4}\ \psi\ =\ \nonumber
\end{equation}
Where the tildes on the $\tilde{H}$ and $\tilde{p}$ defines them as operator. We can use here the series.
\begin{equation}
\sqrt{1 + \tilde{p}^2} ~~=
\end{equation}
\begin{equation}
\mbox{ $1 + \frac{1}{2}\tilde{p}^2 -
\frac{1}{8}\tilde{p}^4 + \frac{1}{16}\tilde{p}^6 - \frac{5}{128}\tilde{p}^8 + \frac{7}{256}\tilde{p}^{10} - \frac{21}{1024}\tilde{p}^{12} +
\frac{33}{2048}\tilde{p}^{14} - \frac{429}{32768}\tilde{p}^{16} + ....$}\nonumber
\end{equation}
This series is found to be represented by a finite expression and we can define $\tilde{H}$ as following.
\begin{equation}
i \hbar \frac{\partial}{\partial t}\, \psi ~~=~~ \tilde{H}\,\psi ~~=\ \sqrt{\mathbf{\tilde{p}}^2 c^2 + m^2c^4}\ \psi\ =\ \nonumber
\end{equation}
\begin{equation}
=\ \pm \ mc^2\left\{ 1 + \sum_{n=1}^\infty
\frac{(-1)^{n+1}\ (2n-2)!}{n!(n-1)!\ 2^{2n-1}} \left(
\frac{\mathbf{\tilde{p}}^2}{m^2c^2}\right)^{n} \right\} \psi
\end{equation}
Or written out as a differential operator:
\begin{equation}
i \hbar \frac{\partial}{\partial t} \, \psi ~~=~~ \tilde{H}\,\psi ~~=
\end{equation}
\begin{equation}
\pm \ mc^2\left\{ 1
- \sum_{n=1}^\infty \frac{ (2n-2)!}{n!(n-1)!\ 2^{2n-1}} \left(
\frac{\hbar^2}{m^2 c^2}\ \right)^{n} \left(
\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}
+\frac{\partial^2}{\partial z^2} \right)^{n}\right\} \psi \nonumber
\end{equation}
Applying this type of operator series on a wave function $\psi$ amounts to a convolution with a Bessel K function. For instance in the one dimensional case one can show, using the Fourier transform, that:
\begin{equation}
\begin{aligned}
& \sqrt{ \partial_x^2 - m^2}^{\,-1} \psi &=~~ &~~\, K_0(mx) ~*~ \psi \\
\end{aligned}
\end{equation}
Where $*$ denotes convolution, From this result then one can derive:
\begin{equation}
\begin{aligned}
& \sqrt{ \partial_x^2 - m^2}~~ \psi &=~~ &\tfrac{m}{x} K_1(mx) ~*~ \psi \\
\end{aligned}
\end{equation}
This convolution means that $\partial\psi/\partial t$ depends on non-local values of $\psi$
Hans