A hypergeometric identity related to Bessel functions
The coefficient of $x^n$ in $_2F_0(\alpha,\beta;z) {}_2F_0(\alpha,\beta;-z)$, multiplied by $n!$, is $$\begin{aligned} \sum_{k=0}^n (-1)^{n-k}\binom nk &(\alpha)_k(\beta)_k (\alpha)_{n-k}(\beta)_{n-k}\hfill\\ &=(-1)^n (\alpha)_n(\beta)_n\, {}_3F_2\left({-n,\atop}\!{\alpha,\atop 1-\alpha -n,}\, {\beta\atop1-\beta-n}\biggm | 1\right). \end{aligned} $$ The right side can be evaluated by Dixon's theorem to give the identity cited by Johannes Trost.
In "Higher Trancendental Functions", Vol. 1, by A. Erdelyi (ed.), on page 86, equation (4) says $$ _2F_0(α,β;z)\ _2F_0(α,β;-z) = \, _4F_1(α,β,\frac{1}{2}(α+β),\frac{1}{2}(α+β+1);α+β;4z^2), $$ from which your formula can be derived by setting $\alpha=−n$, $\beta=n+1$.
Unfortunately Erdelyi, et al., do not give a proof only references, which I reproduce here:
W.N. Bailey (1928), Proc. London Math. Soc. (2) 28, 242-254
Preece, C.T. (1924), Proc. London Math. Soc. (2) 22, 370-380,
but I cannot say more about these papers and their content.
$a(x) = {}_2F_0(-n,n+2; x/2)$ satisfies the differential equation $$ \left( -{n}^{2}-n \right) a \left( x \right) + \left( 2\,x-2 \right) {\frac {\rm d}{{\rm d}x}}a \left( x \right) +{x}^{2}{\frac { {\rm d}^{2}}{{\rm d}{x}^{2}}}a \left( x \right) =0$$ $b(x) = a(-x)$ satisfies $$ \left( -{n}^{2}-n \right) b \left( x \right) + \left( 2\,x+2 \right) {\frac {\rm d}{{\rm d}x}}b \left( x \right) +{x}^{2}{\frac { {\rm d}^{2}}{{\rm d}{x}^{2}}}b \left( x \right)=0$$ $c(x) = {}_3F_0(1/2,-n,n+1; x^2)$ satisfies $$ \left( -4\,{n}^{2}x-4\,nx \right) c \left( x \right) + \left( -4\,{n}^{2}{x}^{2}-4\,n{x}^{2}+6\,{x}^{2}-4 \right) {\frac {\rm d}{{\rm d}x}}c \left( x \right) +6\,{x}^{3}{\frac {{\rm d}^{2}}{ {\rm d}{x}^{2}}}c \left( x \right) +{x}^{4}{\frac {{\rm d}^{3}}{ {\rm d}{x}^{3}}}c \left( x \right)=0$$ Substitute $c(x) = a(x) b(x)$, use the differential equations for $a(x)$ and $b(x)$, and this simplifies to $0=0$. That is, $c(x)$ and $a(x)b(x)$ satisfy the same third order linear differential equation. There is only a one-dimensional family of analytic solutions of this equation ($0$ being an irregular singular point), and the fact that $a(0) b(0) = c(0) = 1$ completes the proof.