Your integral has an elementary closed form, that was correctly stated by Cleo in her answer without proof:
$$S=\int_0^1\left({_2F_1}\!\left(-\tfrac14,\tfrac54;\,1;\,\tfrac{x}2\right)\right)^2dx=\frac{8\sqrt2+4\ln\left(\sqrt2-1\right)}{3\pi}.\tag1$$
Proof: Using DLMF 14.3.6 we can express the hypergeometric function in the integrand as the Legendre function of the $1^{st}$ kind (also known as the Ferrers function of the $1^{st}$ kind) with fractional index:
$${_2F_1}\!\left(-\tfrac14,\tfrac54;\,1;\,\tfrac{x}2\right)=P_{\small1/4}(1-x).\tag2$$
Now the integral can be written as
$$S=\int_0^1\left(P_{\small1/4}(1-x)\right)^2dx=\int_0^1\left(P_{\small1/4}(x)\right)^2dx.\tag3$$
To evaluate it, we use formula 7.113 on page 769 in Gradshteyn & Ryzhyk:
$$\int_0^1P_\nu(x)\,P_\sigma(x)\,dx=\\\frac{\frac{\Gamma\left(\frac12+\frac\nu2\right)\,\Gamma\left(1+\frac\sigma2\right)}{\Gamma\left(\frac12+\frac\sigma2\right)\,\Gamma\left(1+\frac\nu2\right)}\sin\!\left(\frac{\pi\sigma}2\right)\cos\!\left(\frac{\pi\nu}2\right)-\frac{\Gamma\left(\frac12+\frac\sigma2\right)\,\Gamma\left(1+\frac\nu2\right)}{\Gamma\left(\frac12+\frac\nu2\right)\,\Gamma\left(1+\frac\sigma2\right)}\sin\!\left(\frac{\pi\nu}2\right)\cos\!\left(\frac{\pi\sigma}2\right)}{\frac\pi2(\sigma-\nu)(\sigma+\nu+1)}.\tag4$$
Note that in our case $\nu=\sigma=\frac14$, so we cannot use this formula directly because of the term $(\sigma-\nu)$ in the denominator. Instead, we let $\nu=\frac14$ and find the limit for $\sigma\to\frac14$:
$$S=\lim\limits_{\sigma\to{\small1/4}}\int_0^1P_{\small1/4}(x)\,P_\sigma(x)\,dx=\\\lim\limits_{\sigma\to{\small1/4}}\frac{\frac{\Gamma\left(\frac58\right)\,\Gamma\left(1+\frac\sigma2\right)}{\Gamma\left(\frac12+\frac\sigma2\right)\,\Gamma\left(\frac98\right)}\sin\!\left(\frac{\pi\sigma}2\right)\cos\!\left(\frac\pi8\right)-\frac{\Gamma\left(\frac12+\frac\sigma2\right)\,\Gamma\left(\frac98\right)}{\Gamma\left(\frac58\right)\,\Gamma\left(1+\frac\sigma2\right)}\sin\!\left(\frac\pi8\right)\cos\!\left(\frac{\pi\sigma}2\right)}{\frac\pi2(\sigma-\frac14)(\sigma+\frac54)}.\tag5$$
To evaluate the limit, we use l'Hôpital's rule. This gives quite a big expression that I will not copy here. It contains values of the gamma and digamma functions at rational points, that could be simplified to elementary terms using the Gauss digamma theorem and identities given in the MathWorld and in the famous Vidūnas paper, yielding the desired result $(1)$.
Indeed, we can have a more general result:
$$\int_0^1\left({_2F_1}\!\left(-\nu,\nu+1;\,1;\,\tfrac x2\right)\right)^2dx=\int_0^1\left(P_\nu\left(x\right)\right)^2dx=\\\frac{1+\!\left[\psi_0\!\left(1+\frac\nu2\right)-\psi_0\!\left(\frac12+\frac\nu2\right)\right]\frac{\sin(\pi\nu)}\pi}{1+2\nu}.\tag6$$
Consider the hypergeometric equation with parameters $(a,b,c)=\left(\frac16,\frac12,\frac13\right)$, and build from its two canonical solutions near $z=0$ the vector
$$\vec{y}(z)=\left(\begin{array}{c}
y_1 \\ y_2
\end{array}\right)=\left(\begin{array}{c}
_2F_1(a,b;c;z) \\ z^{1-c}{}_2F_1(a-c+1,b-c+1;2-c;z)
\end{array}\right).\tag{1}$$
This is a single-valued vector function on $\mathbb{C}\backslash\{(-\infty,0]\cup[1,\infty)\}$. Its analytic continuation along a closed loop $\gamma$ gives rise to monodromy representation of $\pi_1(\mathbb{C}\backslash\{0,1\})$:
$$ \gamma\mapsto M_{[\gamma]},\qquad y(\gamma z)=M_{[\gamma]}y(z).$$
The monodromy group $G\subset GL(2,\mathbb{C})$ of the hypergeometric equation is generated by two matrices corresponding to simple loops around $0$ and $1$. In the case we are interested in these matrices are explicitly given by
$$M_0=\left(\begin{array}{cc} 1 & 0 \\ 0 & e^{-2\pi i /3}\end{array}\right),\qquad
M_1=C\left(\begin{array}{cc} 1 & 0 \\ 0 & e^{2\pi i /3}\end{array}\right)C^{-1},\tag{2}$$
where the connection matrix $C=\left(\begin{array}{cc} 1 & 2^{\frac43} \\ -2^{\frac83} & 8\end{array}\right)$. If $G$ is finite, then $\vec{y}(z)$ has a finite number of branches, and moreover (Schwarz 1872), is algebraic.
It is not difficult to check that the monodromy group $G$ generated by $M_0$, $M_1$ from (2) is indeed finite. In particular, note that
$$M_0^3=M_1^3=I,\qquad M_1^2=-M_0M_1M_0, $$ $$M_1M_0M_1=-M_0^2,\qquad M_1M_0^2M_1=M_0^2M_1M_0^2.$$
It turns out that $G$ has order $24$ and is isomorphic to the binary tetrahedral group:
$$G\cong 2T=\langle s,t\,|\,(st)^2=s^3=t^3\rangle, $$
where the generators can be identified as $s=M_0M_1M_0M_1M_0$, $t=M_1M_0M_1M_0M_1$.
Corollary: The hypergeometric functions in (1) are algebraic.
Algebraic solutions of the hypergeometric equations are classified by the so-called Schwarz table, and have been studied by many mathematicians, see e.g. the bibliography in this paper by R.Vidunas. Their explicit construction is somewhat involved but relatively straightforward - at least when the corresponding algebraic curve has genus $0$ (the genus can be determined independently from the Riemann-Hurwitz formula).
In our case the task simplifies even more as our parameter values can be obtained from the genus $0$ tetrahedral formula (2.4) of the above mentioned paper by a combination of a linear trasformation (sending $\frac56$ to $\frac43-\frac56=\frac12$) and differentiation (transforming $\frac43$ into $\frac13$). The result is
$$_2F_1\left(\frac16,\frac12;\frac13;-\frac{r(r+2)^3}{(r+1)(1-r)^3}\right)=\frac{\sqrt{1-r^2}}{2r+1}.$$
Corollary: The antiderivative $\displaystyle\int \mathcal{R}\left(x,y(x)\right)dx$, where $y(x)={}_2F_1\left(\frac16,\frac12;\frac13;-x\right)$ and $\mathcal{R}(x,y)$ is rational in both arguments, can be expressed in terms of elementary functions.
Example:
The transformation $r\mapsto x(r)=\frac{r(r+2)^3}{(r+1)(1-r)^3}$ bijectively maps $(0,1)$ to $(0,\infty)$, and therefore the initial integral becomes
\begin{align}\mathcal{I}&=\int_0^1 \left(\frac{\sqrt{1-r^2}}{2r+1}\right)^{12}\left(\frac{r(r+2)^3}{(r+1)(1-r)^3}\right)'dr=\\&=
2\int_0^1 \frac{(1+r)^4(1-r)^2(r+2)^2}{(2r+1)^{10}} dr=\\&=\frac{80\,663}{153\,090}.
\end{align}
Best Answer
Before attacking the integral, I mention something about cubic theta function. The whole solution heavily exploits tools from modular forms. The "footnote" contains more information.
The three cubic theta functions are defined by $$\begin{aligned} a(q) &= \sum_{m,n} q^{m^2+mn+n^2}\\ b(q) &= \sum_{m,n} \zeta_3^{m-n} q^{m^2+mn+n^2}\\ c(q) &= \sum_{m,n} q^{{(m+\frac{1}{3})^2+(m+\frac{1}{3})(n+\frac{1}{3})+(n+\frac{1}{3})^2}} \end{aligned}$$ where $\zeta_3 = e^{2\pi i/3}$, sum is over all $m,n\in \mathbb{Z}$. Then it can be shown$^1$ that $$a(q)^3 = b(q)^3+c(q)^3$$ $$a(q) = \frac{\eta^3(q) + 9 \eta^3(q^9)}{\eta (q^3)}\qquad b(q) = \frac{\eta^3(q)}{\eta(q^3)}\qquad c(q) = 3\frac{\eta^3(q^3)}{\eta(q)}$$ where $\eta(q) = q^{1/24} \prod_{n\geq 1}(1-q^n)$ is the Dedekind eta function.
Define $$K_3(m) = {_2F_1}(\frac{1}{3},\frac{2}{3};1;m) $$ Similar to elliptic integrals, denote $K_3'(m) = K_3(1-m), m' = 1-m$. Then one easily shows (I omit the subscript $3$): $$\frac{d}{dm}(\frac{K'}{K}) = -\frac{\sqrt{3}}{2\pi}\frac{1}{mm'K^2}$$
Moreover, letting $q= \exp(-\frac{2\pi}{\sqrt{3}}\frac{K'(m)}{K(m)})$, the following inversion formula holds$^2$ when $0<m<1$: $$a(q) = K(m)\qquad b(q)=(1-m)^{1/3} K(m)\qquad c(q) = m^{1/3} K(m)$$
Now we tackle the integral, $$I = \frac{1}{3}\int_0^1 {{m^{ - 2/3}}K{{(m)}^2}dm} $$ we make the substitution $q = \exp ( - \frac{{2\pi }}{{\sqrt 3 }}\frac{{K'(m)}}{{K(m)}})$, the above formulas imply $dq = \frac{{q}}{{mm'{K^2}}}dm$, as $m$ increases from $0$ to $1$, $q$ increases from $0$ to $1$. $$I = \frac{1}{3}\int_0^1 {\frac{{b{{(q)}^3}c(q)}}{{mm'{K^2}}}dm} = \frac{1}{3}\int_0^1 {\frac{{b{{(q)}^3}c(q)}}{q}dq} = \int_0^1 {\frac{{\eta {{(q)}^8}}}{q}dq} $$ Next, I will use notation $\eta(q),\eta(\tau)$ interchangably (the common notation in context of modular forms), where $q = e^{2\pi i \tau}$. Make $q=e^{-2\pi x}$, then $I$ becomes $$I = 2\pi \int_0^\infty {\eta {{(ix)}^8}dx} = 2\pi \int_0^\infty {{x^2}\eta {{(ix)}^8}dx} $$ where in last step, I used $\eta(-1/\tau) = \sqrt{-i\tau} \eta(\tau)$. Transform it back to $q$: $$\tag{1} I = \frac{1}{{4{\pi ^2}}}\int_0^1 {\frac{{{{\ln }^2}q}}{q}\eta {{(q)}^8}dq} $$ It can be shown that$^3$: $$\eta {(q)^8} = - \frac{1}{2}\sum\limits_{v \in S} {({v_0} - {v_1})({v_1} - {v_2})({v_0} - {v_2}){q^{{{\left\| v \right\|}^2}/6}}}$$ $$S = \left\{ {v \in {\mathbb{R}^3}|v = ({v_0},{v_1},{v_2}) = (3n,3m + 1,3r - 1),n + m + r = 0,n,m,r\in\mathbb{Z}} \right\}$$ with $\|v\|$ the norm of a vector. Plug this into (1): $$I = \frac{{ - 1}}{{{{(2\pi )}^2}}}{6^3}\sum\limits_{v \in S} {\frac{{({v_0} - {v_1})({v_1} - {v_2})({v_0} - {v_2})}}{{{{\left\| v \right\|}^6}}}} $$ Denote $\rho = e^{\pi i/3}$. Note that $({v_0} - {v_1})({v_1} - {v_2})({v_0} - {v_2}) = 2\Re {({v_0} + \rho {v_1})^3}$ and $${\left\| v \right\|^6} = 8{({v_0} + \rho {v_1})^3}{({v_0} + {\rho ^{ - 1}}{v_1})^3}$$ we obtain $$I = \frac{{ - 27}}{{2{\pi ^2}}}\Re \sum\limits_{v\in S} {\frac{1}{{{{({v_0} + {\rho ^{ - 1}}{v_1})}^3}}}} = - \frac{{27}}{{2{\pi ^2}}}\Re \sum\limits_{(m,n) \in {\mathbb{Z}^2}} {\frac{1}{{{{(3n + {\rho ^{ - 1}}(3m + 1))}^3}}}}$$ The latter can be recognized as an Eisenstein series of level $3$, but to calculate its value, it is best to use Weierstrass elliptic function. Let $\wp_{1,\rho}$ denote this elliptic function with periods $\{1,\rho\}$, then $${\wp _{1,\rho }}'(z) = - 2\sum\limits_{n,m} {\frac{1}{{{{(z + n + m\rho )}^3}}}} $$ gives $$I=\frac{1}{{4{\pi ^2}}}\Re \left[{\wp _{1,\rho }}'(\frac{{{\rho ^{ - 1}}}}{3})\right] = \frac{{{\omega ^3}}}{{4{\pi ^2}}}\Re\left[ {\wp _{\omega ,\omega \rho }}'(\frac{{{\omega\rho ^{ - 1}}}}{3})\right]$$ where $\omega = \Gamma(1/3)^3/(2\pi)$, then it is well-known that modular invariants associated to periods $\{\omega,\omega\rho\}$ are $g_2 = 0, g_3 = 1$. Therefore ${\wp _{\omega ,\omega \rho }}'(\frac{{\omega {\rho ^{ - 1}}}}{3})$ is the $y$-coordinate of a $3$-torsion of the elliptic curve $y^2 = 4x^3 - g_2 x - g_3 = 4x^3 -1$, which can be readily calculated to be $\sqrt{3}$. Finally we finish the calculation:$I = \omega^3\sqrt 3/(4\pi^2)$.
$^1$: Proof outline: $a(q^3),b(q^3),c(q^3)$ are modular forms of weight $1$ and level $27$, therefore it suffices to verify their $q$-expansions to certain power of $q$. A self-contained approach can be found in the 1994 paper Cubic Analogues of the Jacobian Theta Function.
$^2$: Proof outline: $f=c^3(\tau)/a^3(\tau)$ is modular function of $\Gamma_0(3)$, by a fact in modular forms, $b(\tau)$ satisfies a 2nd order ODE in terms of $f$, its coefficients are rational functions of $f$ since modular curve $X(3)$ has genus $0$. Therefore, in certain region of $\mathbb{H}$, $b(\tau) = (1-f)^{1/3} K_3(f)$, we could replace $\tau$ by $\gamma\tau$ for $\gamma\in \Gamma_0(3)$, modularity of $b$ allows us to isolate the $\tau$. But doing this replacement might change it into another linear independent solution of the ODE, which explains why $K'/K$ arises. The details are more delicate.
$^3$: The exponent $8$ is special here, which is the dimension of semisimple Lie algebra $A_2$. There is a corresponding formula for $\eta(q)^d$ each semisimple Lie algebra with dimension $d$. See Affine Root Systems and Dedekind's eta-Function by I.G. Macdonald.