Encontré una demostración elemental usando fracciones parciales. Esto está motivado por la observación de Nemo de que series muy similares aparecen en el trabajo de Schlosser (Sección 7). La inversión de matriz utilizada por Schlosser está estrechamente relacionada con las fracciones parciales, ver por ejemplo este artículo.
Fijemos $\lambda$ y sean $t_1$ y $t_2$ variables que satisfacen $$t_1+t_2-\frac{t_1t_2}{\lambda}=1-\frac{1}{4\lambda}. $$ Por hechos elementales sobre expansiones en fracciones parciales, podemos escribir $$ \label{pf}\frac{(t_1+t_2+1)_n}{(t_1)_{n+1}(t_2)_{n+1}}=\sum_{k=0}^n A_k\left(\frac{1}{t_1+k}+\frac 1{t_2+k}-\frac 1{\lambda+k}\right), \tag 1 $$ donde $$ A_k= \frac{(t_1+k)(t_1+t_2+1)_n}{(t_1)_{n+1}(t_2)_{n+1}}\Bigg|_{t_1=-k}. $$ El último término en \eqref{pf} hace que cada término se anule en el límite $t_1\rightarrow\infty$, $t_2\rightarrow \lambda$ (y viceversa). Escribiendo $$ a_k=t_1+t_2\Big|_{t_1=-k}=1-\frac{(2k+1)^2}{4(\lambda+k)}, $$ esto se puede simplificar a $$ A_k=\frac{(a_k+1)_{k-1}(-1)^k}{k!(n-k)!(a_k+n+1)_{k}}. $$ Especializando $t_1=t_2=1/2$ en \eqref{pf} obtenemos $$\frac{(2)_n}{(1/2)_{n+1}^2}=\sum_{k=0}^n \frac{(a_k+1)_{k-1}(-1)^k}{k!(n-k)!(a_k+n+1)_{k}}\left(\frac{4}{2k+1}-\frac 1{\lambda+k}\right). $$ Ahora multiplicamos ambos lados por $n!$ y dejamos que $n\rightarrow\infty$. El lado izquierdo se puede escribir $$\frac{\Gamma(n+2)\Gamma(n+1)}{\Gamma(n+3/2)^2}\cdot\frac{\Gamma(1/2)^2}{\Gamma(2)}\rightarrow \frac{\Gamma(1/2)^2}{\Gamma(2)}=\pi.$$ En el lado derecho, tenemos el factor $$\frac{n!}{(n-k)!(a_k+n+1)_{k}}=(-1)^{k}\frac{(-n)_{k}}{(a_k+n+1)_{k}}\rightarrow 1. $$ Esto nos da (no debería ser difícil justificar tomar el límite) \begin{align*}\pi&=\sum_{k=0}^\infty\frac{(-1)^k}{k!}\left(\frac{4}{2k+1}-\frac 1{\lambda+k}\right)\left(2-\frac{(2k+1)^2}{4(\lambda+k)}\right)_{k-1}\\ &=\sum_{k=0}^\infty\frac{1}{k!}\left(\frac 1{\lambda+k}-\frac{4}{2k+1}\right)\left(\frac{(2k+1)^2}{4(\lambda+k)}-k\right)_{k-1} , \end{align*} que es la identidad deseada.
Más generalmente, se podrían tomar variables que satisfagan $$t_1+t_2-\frac{t_1t_2}{\lambda}=x+y-\frac{xy}{\lambda}, $$ tomar la expansión en fracciones parciales de $$\frac{(t_1+t_2-p)_n}{(t_1)_{n+1}(t_2)_{n+1}}, $$ especializar $t_1=x$, $t_2=y$ y finalmente dejar que $n\rightarrow\infty$. Comprobé que esto da la ecuación (4) en el artículo de Sana y Sinha (donde $x=\alpha-s_1$, $y=\alpha-s_2$). Presumiblemente otras identidades en su artículo siguen por variaciones adicionales del mismo argumento.