Integral's Closed-form expression in terms of hypergeometric function I want to solve the following integral:

$$I = 2\left[\int_{0}^{1}\dfrac{y^m}{(1 - ay)^{m + 1}\sqrt{1 - y^2}}\mathrm{d}y+\int_{0}^{1}\dfrac{y^m}{(1 + ay)^{m + 1}\sqrt{1 - y^2}}\mathrm{d}y\right]$$

where $0 \leq a \leq 1$ and $m$ is an arbitrary positive constant.
I have tried to expand each term of the above integral in series form and solve it. I have got answer in an infinite series form, which is not a closed-form expression. 
Infinite series expression (Derivation is given below):
$$
\boxed{
I = 4\sum_{t = 0}^{\infty}\sum_{\underset{u \neq \text{Odd}}{u=0}}^{t}\frac{a^u(m+1)_{u}}{u!}\frac{(1/2)_{t-u}}{(t-u)!}\frac{1}{(m+1)+2t-u}.}
$$
However, Maple gives a beautiful closed-form expression in terms of hypergeometric function, which is as follows:
$$I = \frac{2^{m+1}\left[\sqrt{\pi}\Gamma\left(\frac{m + 1}{2}\right)^2 {_2F_1\left(\left[\frac{m + 1}{2},\frac{m+ 1}{2}\right], \left[\frac{1}{2}\right], a^2\right)}\right]}{\sqrt{\pi}\Gamma(m + 1)}$$
where ${_2F_1\left(\left[a, b\right], \left[w\right], z \right)}$ is a hypergeometric function.
How can I obtain the above expression given by Maple? 
Derivation of infinite series expression: Following is the derivation for the infinite series expression that I could figure out.
Let the first integral of $I$  be
$$I_1 = \int_{0}^{1}\dfrac{y^m}{(1 - ay)^{m + 1}\sqrt{1 - y^2}}\mathrm{d}y = \int_{0}^{1} y^m(1 - ay)^{-(m + 1)}(1 - y^2)^{-1/2}\mathrm{d}y.$$ 
Then, we make use of the following relation
$$ (1 - y)^{-m} = \sum_{t=0}^{\infty}\frac{(m)_{t}}{t!}y^{t}$$
where $(m)_{t}$ denotes the Pochhammer symbol and is defined as $(m)_{t} = (t)\times(m+1)\times\dotsc\times(m+t-1) = \frac{\Gamma(m+t)}{\Gamma(m)}$.
Then, we have
$$ (1 - ay)^{-(m+1)} = \sum_{t=0}^{\infty}\frac{(m+1)_{t}}{t!}(ay)^{t}.$$
On the similar lines, one can write 
$$(1 - y^2)^{-1/2} = \sum_{u=0}^{\infty}\frac{(1/2)_{u}}{u!}y^{2u}.$$
Thus, $I_1$ can be written as
$$I_1 =  \int_{0}^{1} y^m\left(\sum_{t=0}^{\infty}\frac{a^t(m+1)_{t}}{t!}y^{t}\right)\left(\sum_{u=0}^{\infty}\frac{(1/2)_{u}}{u!}y^{2u}\right)\mathrm{d}y.$$
Simplifying further,
$$I_1 = \int_{0}^{1} \left(\sum_{t=0}^{\infty}\frac{a^t(m+1)_{t}}{t!}y^{t+m}\right)\left(\sum_{u=0}^{\infty}\frac{(1/2)_{u}}{u!}y^{2u}\right)\mathrm{d}y.$$
We make use of the following relation
$$
\left(\sum_{p=0}^{\infty}z_p\right)\left(\sum_{q=0}^{\infty}g_q\right) = \sum_{t = 0}^{\infty}\sum_{u=0}^{t}z_u g_{u-t}.
$$
Then, we can write
$$
I_1 = \int_{0}^{1} \sum_{t = 0}^{\infty}\sum_{u = 0}^{t}\frac{a^uy^{u+m}(m+1)_{u}}{u!}\frac{y^{2(t-u)}(1/2)_{t-u}}{(t-u)!}\mathrm{d}y.$$
Then,
$$
I_1 = \int_{0}^{1} \sum_{t = 0}^{\infty}\sum_{u = 0}^{t}\frac{a^u(m+1)_{u}}{u!}\frac{(1/2)_{t-u}}{(t-u)!}y^{(m+2t-u)}\mathrm{d}y.
$$
After integration, we can write
$$
I_1 = \sum_{t = 0}^{\infty}\sum_{u = 0}^{t}\frac{a^u(m+1)_{u}}{u!}\frac{(1/2)_{t-u}}{(t-u)!}\frac{1}{(m+1)+2t-u}.
$$
Similarly, 
$$
I_2 = \int_{0}^{1}\dfrac{y^m}{(1 + ay)^{m + 1}\sqrt{1 - y^2}}\mathrm{d}y = \sum_{t = 0}^{\infty}\sum_{u = 0}^{t}\frac{(-1)^u a^u(m+1)_{u}}{u!}\frac{(1/2)_{t-u}}{(t-u)!}\frac{1}{(m+1)+2t-u}.
$$
Thus, we can write $I$ as
$$
\boxed{
I = 4\sum_{t = 0}^{\infty}\sum_{\underset{u \neq \text{Odd}}{u=0}}^{t}\frac{a^u(m+1)_{u}}{u!}\frac{(1/2)_{t-u}}{(t-u)!}\frac{1}{(m+1)+2t-u}.}
$$
 A: Expand only the factors $(1\pm ay)^{-m-1}$ in powers of $a$ using
$$(1\mp ay)^{-m-1}=\sum_{k=0}^{\infty}\frac{(m+1)_k}{k!}\left(\pm ay\right)^k,$$
where $(z)_k:=\displaystyle \frac{\Gamma(z+k)}{\Gamma(z)}$ denotes the Pochhammer's symbol. 
The odd powers clearly cancel out and even ones are equal for the two summands, hence we can write
$$I=4\sum_{k=0}^{\infty}\frac{(m+1)_{2k}}{(2k)!} a^{2k}\int_0^{1}\frac{y^{2k+m}dy}{\sqrt{1-y^2}}.$$
The change of variables $y^2=x$ reduces the remaining integral to Euler beta function:
$$\int_0^{1}\frac{y^{2k+m}dy}{\sqrt{1-y^2}}=\frac12 \int_0^{1}\frac{x^{k+\frac{m-1}{2}}dx}{\sqrt{1-x}}=\frac12B\left(k+\frac{m+1}{2},\frac12\right)=\frac{\sqrt\pi}{2}\frac{\Gamma\left(k+\frac{m+1}{2}\right)}{
\Gamma\left(k+1+\frac{m}{2}\right)},$$
so that the initial integral becomes
$$I=2\sqrt\pi \sum_{k=0}^{\infty}\frac{(m+1)_{2k}}{(2k)!}\frac{\Gamma\left(k+\frac{m+1}{2}\right)}{
\Gamma\left(k+1+\frac{m}{2}\right)} a^{2k}.\tag{1}$$
Let us now rewrite the expression for the coefficient of this series using the gamma function duplication formula $\Gamma(2z)=\displaystyle\frac{2^{2z-1}}{\sqrt\pi}\Gamma\left(z\right)\Gamma\left(z+\frac12\right)$. We get
\begin{align*}
\frac{(m+1)_{2k}}{(2k)!}\frac{\Gamma\left(k+\frac{m+1}{2}\right)}{
\Gamma\left(k+1+\frac{m}{2}\right)}&=\frac{\color{red}{\Gamma\left(2k+m+1\right)}\Gamma\left(k+\frac{m+1}{2}\right)}{\color{red}{\Gamma\left(2k+1\right)}\Gamma\left(m+1\right)
\Gamma\left(k+1+\frac{m}{2}\right)}=\\
&=\color{red}{2^m\frac{\Gamma\left(k+\frac{m+1}2\right)\Gamma\left(k+1+\frac{m}2\right)}{\Gamma\left(k+\frac12\right)\Gamma\left(k+1\right)}}\frac{\Gamma\left(k+\frac{m+1}{2}\right)}{\Gamma\left(m+1\right)
\Gamma\left(k+1+\frac{m}{2}\right)}=\\
&=\frac{2^m}{\Gamma(m+1)}\frac{\Gamma^2\left(k+\frac{m+1}{2}\right)}{k!\,\Gamma\left(k+\frac12\right)}=\\
&=\frac{2^m\Gamma^2\left(\frac{m+1}{2}\right)}{\sqrt\pi \,\Gamma(m+1)}\cdot
\frac{\left(\frac{m+1}{2}\right)_k\left(\frac{m+1}{2}\right)_k}{k!\,\left(\frac{1}{2}\right)_k}.
\end{align*}
This immediately reduces the series (1) to the defining series representation of the hypergeometric function:
$$ 
\boxed{\quad\displaystyle I=2^{m+1}\frac{\Gamma^2\left(\frac{m+1}{2}\right)}{\,\Gamma(m+1)}\,_2F_1\left[\frac{m+1}{2},\frac{m+1}{2};\frac12;a^2\right]\quad}$$
