We first use the differentiation formula for the generalized hypergeometric function
\begin{equation}
\frac{a_1a_2\dots a_{p}}{b_1b_2\dots b_q}{}_pF_q\left(\left.\begin{array}{c} c_1,c_2,\dots ,c_p\\ d_1,d_2,\dots ,d_q \end{array}\right| z\right)=\frac{d}{dz}{}_pF_q\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_p\\ b_1,b_2,\dots ,b_q \end{array}\right| z\right)
\end{equation}
Then, the LHS of the proposed identity can be written as
\begin{equation}
_pF_q\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_p\\ b_1,b_2,\dots ,b_q \end{array}\right| z\right)+z\,\frac{a_1a_2\dots a_{p-1}}{b_1b_2\dots b_q}{}_pF_q\left(\left.\begin{array}{c} c_1,c_2,\dots ,c_p\\ d_1,d_2,\dots ,d_q \end{array}\right| z\right)=\left( 1+\frac{z}{a_p}\frac{d}{dz} \right){} _pF_q\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_p\\ b_1,b_2,\dots ,b_q \end{array}\right| z\right)\tag{1}\label{eq1}
\end{equation}
To differentiate the hypergeometric function, we use the Euler's integral transform
\begin{align}
& _pF_q\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_p\\ b_1,b_2,\dots ,b_q \end{array}\right| z\right)\\
&=\frac{\Gamma(b_q)}{\Gamma(a_p)\Gamma(b_q-b_p)} \int_0^1t^{a_p-1}\left( 1-t \right)^{b_q-a_p-1}{}_{p-1}F_{q-1}\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_{p-1}\\ b_1,b_2,\dots ,b_{q-1} \end{array}\right| t\right)\,dt
\end{align}
Here $b_q=a_p+1$, then
\begin{align}
_pF_q\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_p\\ b_1,b_2,\dots ,b_q \end{array}\right| z\right)&=
a_p \int_0^1t^{a_p-1}{}_{p-1}F_{q-1}\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_{p-1}\\ b_1,b_2,\dots ,b_{q-1} \end{array}\right| zt\right)\,dt\\
&=\frac{a_p}{z^{a_p}} \int_0^zu^{a_p-1}{}_{p-1}F_{q-1}\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_{p-1}\\ b_1,b_2,\dots ,b_{q-1} \end{array}\right| u\right)\,du
\end{align}
Then
\begin{align}
\frac{d}{dz}&\,{} _pF_q\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_p\\ b_1,b_2,\dots ,b_q \end{array}\right| z\right)\\
&=\frac{a_p}{z}\,{}_{p-1}F_{q-1}\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_{p-1}\\ b_1,b_2,\dots ,b_{q-1} \end{array}\right| z\right)-\frac{a_p}{z} \,{}_pF_q\left(\left.\begin{array}{c} a_1,a_2,\dots ,a_p\\ b_1,b_2,\dots ,b_q \end{array}\right| z\right)
\end{align}
Plugging this expression in eq. \eqref{eq1} we find theRHS of the proposed identity.
Best Answer
My question is related to this question by Vladimir. Because it is already proved that $${\large\int}_{-1}^1\frac{dx}{\sqrt[3]{9+4\sqrt5\,x}\ \left(1-x^2\right)^{\small2/3}} = \frac{3^{\small3/2}}{2^{\small4/3}\,5^{\small5/6}\,\pi }\Gamma^3\!\!\left(\tfrac13\right),$$
we have the answer for my question too, since according to Maple
$${\large\int}_{-1}^1\frac{dx}{\sqrt[3]{9+4\sqrt5\,x}\ \left(1-x^2\right)^{\small2/3}} = \frac{2}{9} \frac{\sqrt[3]{4}\,\sqrt[3]{3}\,\pi^2\,{_2F_1}\left(\begin{array}c\tfrac16,\tfrac23\\\tfrac56\end{array}\middle|\,\frac{80}{81}\right)}{\Gamma^3 \left( \frac{2}{3} \right)}.$$
The second part of the question $-$ the part under the line $-$ is still open.