In my research, I came across the following formula for Apéry numbers:
$$ a_n = \sum_{i=0}^n \sum_{j=0}^n {n \choose i}^2 {n \choose j}^2{i+j \choose i}.$$
This formula does not seem to appear in the OEIS entry for Apéry numbers, despite the entry containing several other expansions as binomial coefficients sums. Does this formula appear anywhere in the literature on Apéry numbers?
It seems a bit restrictive to look for a reference only "in the literature on Apéry numbers". I found a reference in the much larger literature on hypergeometric functions. In terms of Pochhammer symbols $$(a)_k=a(a+1)\dotsm(a+k-1),$$ the identity $a_n=b_n$ takes the form $$\sum_{i,j=0}^n\frac{(1)_{i+j}(-n)_i(-n)_i(-n)_j(-n)_j}{(1)_i(1)_ii!(1)_j(1)_jj!}=\sum_{k=0}^n\frac{(-n)_k(-n)_k(n+1)_k(n+1)_k}{(1)_k(1)_k(1)_kk!}$$ or in hypergeometric notation $$F^{1:2}_{0:2}\left(\begin{matrix}1&:&-n,-n&;&-n,-n\\ &:&1,1&;&1,1\end{matrix};1,1\right)={}_4F_3\left(\begin{matrix}-n,-n,n+1,n+1\\1,1,1\end{matrix};1\right).$$ Reversing the summation on the right-hand side gives the alternative version $$F^{1:2}_{0:2}\left(\begin{matrix}1&:&-n,-n&;&-n,-n\\ &:&1,1&;&1,1\end{matrix};1,1\right)=\binom{2n}n^2{}_4F_3\left(\begin{matrix}-n,-n,-n,-n\\1,-2n,-2n\end{matrix};1\right).$$ This can be recognized as the special case $a=a'=b'=-n$, $c=c'=e=1$ of the transformation formula $$F^{1:2}_{0:2}\left(\begin{matrix}e&:&a,-n&;&a',b'\\ &:&c,e&;&c',e\end{matrix};1,1\right)=\frac{(c-a)_n\Gamma(c')\Gamma(c'-a'-b')}{(c)_n\Gamma(c'-a')\Gamma(c'-b')}{}_4F_3\left(\begin{matrix}-n,a,a',b'\\e,1+a-c-n,1+a'-c'+b'\end{matrix};1\right)$$ given as equation (22) in J. Van der Jeugt, S. N. Pitre and K. Srinivasa Rao, Transformation and summation formulas for double hypergeometric series, J. Comput. Appl. Math. 83 (1997), 185–193.
Alternatively, one can reverse the summation on both sides, so that the left-hand side becomes a series of the form $F^{0:3}_{1:1}$. Then, the identity $a_n=b_n$ follows from equation (5) in P. W. Karlsson, Two hypergeometric summation formulae related to 9-j coefficients, J. Phys. A 27 (1994), 6943–6945.
Both papers mentioned are motivated by $9j$-symbols of quantum mechanics. This is not a coincidence. As discussed in R. Askey and J. A. Wilson, A recurrence relation generalizing those of Apéry, J. Austral. Math. Soc. Ser. A 36 (1984), 267–278, the recursions used by Apéry in his proof of the irrationality of $\zeta(3)$ are special cases of the recursions for Wilson polynomials, which appear as $6j$-symbols of the group $\mathrm{SU}(1,1)$. That is, Apéry numbers are $6j$-symbols. The more general $9j$-symbols are given by hypergeometric triple sums, which reduce to double sums like $F^{1:2}_{0:2}$ in special cases. Transformations like the one you ask about appear from even more special situations when the $9j$-symbol in fact equals a $6j$-symbol.
I'm not sure about references, but we can show that this formula and the title formula in the OEIS entry give the same numbers.
The formula in the OEIS reads: $$b_n:=\sum_{k=0}^n \binom{n}k^2 \binom{n+k}k^2,$$ and below we show that $a_n=b_n$ for all $n$.
First, we note that for a fixed $n$ the generating function for $\binom{n}k^2$ expresses in terms of Legendre polynomials as $$\sum_{k\geq 0} \binom{n}k^2 t^k = (1-t)^n P_n\left(\frac{1+t}{1-t}\right).$$ Then, since $\binom{i+j}{i} = [y^i]\ (1-y)^{-j-1}$, we have $$\sum_{i=0}^n \binom{n}i^2 \binom{i+j}{i} = [y^n]\ (1-y)^{n-j-1} P_n\left(\frac{1+y}{1-y}\right)$$ and thus \begin{split} a_n &= [y^n]\ \sum_{j=0}^n (1-y)^{n-j-1} P_n\left(\frac{1+y}{1-y}\right)\cdot [t^j]\ (1-t)^n P_n\left(\frac{1+t}{1-t}\right)\\ &=[y^n]\ (1-y)^{n-1} P_n \left(\frac{1+y}{1-y}\right)\cdot \left(\frac{y}{1-y}\right)^n P_n \left(\frac{2-y}{y}\right) \\ &= [y^0]\ (1-y)^{-1} P_n \left(\frac{1+y}{1-y}\right) P_n \left(\frac{2-y}{y}\right). \end{split} Finally, using the formulas for Legendre polynomials: $$ P_n(\frac{1+y}{1-y}) = \sum_{k=0}^n \binom{n}k \binom{n+k}k \left(\frac{y}{1-y}\right)^k$$ and $$ P_n(\frac{2-y}{y}) = \sum_{j=0}^n \binom{n}j \binom{n+j}j \left(\frac{1-y}{y}\right)^j,$$ we conclude that \begin{split} a_n &= [y^0]\ (1-y)^{-1} P_n \left(\frac{1+y}{1-y}\right) P_n \left(\frac{2-y}{y}\right) \\ &= [y^0] \sum_{k=0}^n \binom{n}k \binom{n+k}k \sum_{j=0}^n \binom{n}j \binom{n+j}j y^{k-j} (1-y)^{j-k-1} \\ &=\sum_{k=0}^n \binom{n}k^2 \binom{n+k}k^2\\ &=b_n, \end{split} where we noticed that $[y^0]\ y^{k-j} (1-y)^{j-k-1}$ is nonzero (equal 1) only if $j=k$.
While there are already two good answers, I have found a more elementary combinatorial proof of the identity $$\sum_{i=0}^n \sum_{j=0}^n {n \choose i}^2 {n \choose j}^2{i+j \choose i} = \sum_{k=0}^n {n\choose k}^2{n+k \choose k}^2 \qquad (1)$$ that I think is worth sharing here.
First rewrite $$ \sum_{i=0}^n \sum_{j=0}^n {n \choose i}^2 {n \choose j}^2{i+j \choose i} = \sum_{i=0}^n \sum_{j=0}^n \sum_{\ell=0}^{\min(i,j)} {n \choose i}^2 {n \choose j}^2{i \choose \ell}{j \choose \ell}=\sum_{\ell=0}^{n} \left(\sum_{i=\ell}^n {n \choose i}^2{i \choose \ell}\right)^2.$$ By the Vandermonde identity, we have $$\sum_{i=\ell}^n {n\choose i}^2{i \choose \ell} = \sum_{i=\ell}^n {n\choose i}{n-\ell \choose n-i}{n\choose \ell}={2n-\ell \choose n}{n\choose \ell},$$ so we get $$\sum_{i=0}^n \sum_{j=0}^n {n \choose i}^2 {n \choose j}^2{i+j \choose i}=\sum_{\ell=0}^n {n\choose \ell}^2{2n-\ell \choose n}^2.$$ Letting $k=n-\ell$ yields the desired identity $(1)$.
Unpacking the binomial coefficient manipulations used above, we can also rephrase this as a double counting proof as follows:
Let $A, B, C, D$ be four disjoint finite sets of size $n$. We will count in two ways the number of tuples $(X,Y,Z)$ where
satisfying
On the one hand, the number of ways to choose $X,Y$ such that $\# (X\cap A)=n-i$ and $\#(Y\cap B) = j$ is ${n \choose i}^2 {n \choose j}^2$. Then $Z \subseteq A\sqcup B$ is a subset of size $n$ and has to be sandwiched between $X\cap A$ (which has size $n-i$) and $A\sqcup (B\cap Y)$ (which has size $n+j$). The number of ways to choose such a $Z$ is ${i+j \choose i}$. Therefore the number of tuples $(X,Y,Z)$ is given by the left-hand side of $(1)$.
On the other hand, there are ${n \choose k}^2$ ways to pick $Z$ such that $\#(Z\cap A)=k$. For every such choice of $Z$, there are ${n \choose n+k}$ ways to choose $X$ satisfying condition 1. and ${k \choose n+k}$ ways to choose $Y$ satisfying condition 2. Therefore the number of tuples $(X,Y,Z)$ is given by the right-hand side of $(1)$.
