1 st proof. By elementary algebra
$$
\binom{n}{x}^2\binom{2n+x}{x}=\frac{(n!)^2}{(2n)!}\frac{\prod_{k=1}^{2n}(x+k)}{\prod_{k=1}^{n}(x-k)^2}\frac{\sin^2 \pi x}{\pi^2 x^2}.
$$
Now one can apply residue theorem to the function $\frac{(n!)^2}{(2n)!}\frac{\prod_{k=1}^{2n}(z+k)}{\prod_{k=1}^{n}(z-k)^2}\frac{1-e^{2\pi i z}}{2\pi^2 z^2}$ which has only simple poles at $z=m,~0\le m\le n$ and asymptotics $O(1/z^2),~z\to\infty$ in the upper half plane. The contour of integration is the interval $(-R,R)$ with small semicircles near the poles, and closed by a semicircle of radius $R$ in the upper half plane.
\begin{align}
\int_{-\infty}^{\infty}\binom{n}{x}^2\binom{2n+x}{x} dx&=\text{Re}\left\{\int_{-\infty}^{\infty}\frac{(n!)^2}{(2n)!}\frac{\prod_{k=1}^{2n}(z+k)}{\prod_{k=1}^{n}(z-k)^2}\frac{1-e^{2\pi i z}}{2\pi^2 z^2}dz\right\}\\&=\text{Re}\left\{\pi i \frac{(n!)^2}{(2n)!}\sum_{m=0}^n{\text{res}}_{z=m}\frac{\prod_{k=1}^{2n}(z+k)}{\prod_{k=1}^{n}(z-k)^2}\frac{1-e^{2\pi i z}}{2\pi^2 z^2}\right\}\\
&=\text{Re}\left\{-2\pi i\frac{i}{2\pi}\frac{(n!)^2}{(2n)!}\left(\frac{(2n)!}{(n!)^2}+\sum_{m=1}^n\frac{\prod_{k=1}^{2n}(m+k)}{(m!(n-m)!)^2}\right)\right\}\\&=\frac{(n!)^2}{(2n)!}\sum_{m=0}^n\frac{\prod_{k=1}^{2n}(m+k)}{(m!(n-m)!)^2}\\&=\sum_{m=0}^{n}\binom{n}{m}^2\binom{2n+m}{m}=\binom{2n}{n}^2
\end{align}
2 nd Proof. This approach gives the following generalization.
If $m,n\in\mathbb{N},~0<\beta<1$
$$
f(x)=\binom{m}{\beta x}\binom{n}{\beta x}\binom{m+n+\beta x}{m+n}
$$
then for arbitrary $s\in\mathbb{C}$
$$
\beta\sum_{k=-\infty}^\infty f(s+k)=\beta\int_{-\infty}^\infty f(x)dx=\binom{m+n}{n}^2.
$$
An outline of the proof is given below. First by observing that
$$
f(x)=\frac{m!n!}{(m+n)!}\frac{\prod_{k=1}^{m+n}(\beta x+k)}{\prod_{k=1}^{m}(\beta x-k)\prod_{k=1}^{n}(\beta x-k)}\frac{\sin^2\pi\beta x}{\pi^2\beta^2 x^2},
$$one can see that $f(x)$ does not have any poles and $f(z)\sim e^{2\pi\beta |\text{Im}z|}/|z|^2$, $z\to\infty$.
Now consider the Fourier transform
$$
\int_{-\infty}^\infty f(x)e^{iax}dx.
$$
Contour integration shows that if $|a|>2\pi\beta$ this integral is $0$. This means that $f(x)$ is band limited, i.e. its Fourier spectrum is limited to the band $|\omega|<2\pi\beta$. Now Poisson summation gives
$$
\sum_{k=-\infty}^\infty f(s+k)=\int_{-\infty}^\infty f(x)dx,\quad 2\pi\beta<2\pi.
$$
This means that the sum on the left does not depend on $s$ and $\beta$. So one can put $s=0$, $\beta=1$ and use the result to compute this sum.$\Box$
