What's wrong with the square root version of the Klein-Gordon equation? The Wikipedia article has a derivation of the Klein-Gordon equation.  It gets to this step:
$$\sqrt{\textbf{p}^2 c^2 + m^2 c^4} = E$$
and inserts the QM operators to get
$$\left( \sqrt{ (-i \hbar \nabla)^2 c^2 + m^2 c^4 } \right) \psi = i\hbar\frac{\partial}{\partial t} \psi$$
The article then says

This, however, is a cumbersome expression to work with because the differential operator cannot be evaluated while under the square root sign. In addition, this equation, as it stands, is nonlocal.

To fix this, the first equation is squared instead to get
$$\textbf{p}^2 c^2 + m^2 c^4 = E^2$$
after which the QM operators are inserted and the expression is simplified to get
$$-\hbar^2 c^2 \nabla^2 \psi + m^2 c^4 \psi = -\hbar^2 \frac{\partial^2}{\partial t^2} \psi$$
A couple things I don't understand.  First, are the solutions to this differential equation not exactly the same as the solutions to the first differential equation?  Both sides of the starting equation were squared, so it seems to me that regardless of the particular form of the resulting differential equation, both of them should have the exact same set of solutions.
Secondly, why is the first differential equation cumbersome to work with?  It seems like it would in fact be easier to work with, since the operator under the square root could be expanded in terms of a Taylor series and then you have an equation that is first order in time.
And finally, can someone elaborate on what nonlocal means? The linked article on the Wikipedia page didn't entirely help me understand it.
 A: 
Secondly, why is the first differential equation cumbersome to work with? It seems like it would in fact be easier to work with, since the operator under the square root could be expanded in terms of a Taylor series and then you have an equation that is first order in time.

Well, the Taylor series for operator expressions only really make sense if they converge everywhere (e.g., $\exp(\partial_{x})$ makes sense as a series expression)...modulo technical details.
The squareroot doesn't have a nice series expression for any operator...it works for normal operators.
So what happens with this squareroot version of the KG equation, we just take the Fourier transform of the expression
$$ \int (k^{2}-m^{2})\,\hat{f}(k)e^{ikx}\,\mathrm{d}^{4}k=0. $$
Then observe we have the operator be $(k^{2}-m^{2})$. So, hey, presto, take its squareroot! We get
$$ \int \sqrt{k^{2}-m^{2}}\;\hat{f}(k)e^{ikx}\,\mathrm{d}^{4}k=0. $$
Then...well, then it's a pain to work with. Why? Because all our lovely tools from linear algebra don't really work too well. My next tool, profanity, doesn't produce much results either :\
Addendum: I thought I ought to add some links on this, because there are people researching it. (This method I sketched describes treating the squareroot of the Klein-Gordon Equation using pseudodifferential operators)


*

*Claus Lämmerzahl, "The pseudodifferential operator square root of the Klein–Gordon equation". J. Math. Phys. 34 9 (1993), 3918-3932, doi:10.1063/1.530015

*J. Sucher, "Relativistic Invariance and the Square‐Root Klein‐Gordon Equation". J. Math. Phys. 4 17 (1963); doi:10.1063/1.1703882
I'm sure from there, you can follow the references to where-ever you want.

And finally, can someone elaborate on what nonlocal means? The linked article on the Wikipedia page didn't entirely help me understand it.

As I understand it (and someone will probably correct me if I am wrong), generically, it means the field at one point depends on its value at other spatially separated points. It borks up our intuitive understanding of cause and effect.
If we have infinitely many derivatives, we get this problem. Why?
Well, consider a special case: the Taylor expansion. We have
$$ f(x+h) = f(x) + hf'(x) +\cdots = \exp(h\partial_{x})f(x) $$
where $\exp(h\partial_{x})=1 + h\partial_{x} + \cdots$ is an expression involving infinitely many derivatives. We then get a relation between values at two distinct points ($x$ and $x+h$).
More generally, we could consider any operator involving infinitely many derivatives, not just $\exp(h\partial/\partial x)$.
A: On the right-hand side we have a square root of an operator.  It is possible to take the square root of an operator (ie. the square root of a matrix) and there is a body of theory in linear algebra and spectral theory related to this possibility but the question is how to interpret this physically, since a matrix has multiple square roots which are themselves matrices. 
One possible interpretation is to make a Taylor series expansion as you rightly say, but we then get a Hamiltonian with derivatives of arbitrarily high order.  The two standard approaches are obviously to square both sides and obtain the Klein-Gordon equation or just to propose a Hamiltonian which is linear in the momentum and equal to the square of the relativistic energy-momentum relation: this leads to the Dirac equation.  If you take the latter approach, a solution of the equation is not just a function only and has to have four components.
