The integral of an elliptic integral: $\int_{0}^{1}\frac{x\mathbf{K}^2\left ( x \right )}{\sqrt{1-x^{2}}}\mathrm{d}x$

This result and a sketch for its proof can be found in Moments of products of elliptic integrals by J.G. Wan. It uses a theorem by Zudilin which gives a triple integral representation for some very-well poised $\,_7F_6$ hypergeometric series : \begin{align} &\iiint_{\left[ 0,1\right]^3}\frac{x^{h_2-1}y^{h_3-1}z^{h_4-1}(1-x)^{h_0-h_2-h_3 }(1-y)^{h_0-h_3-h_4}(1-z)^{h_0-h_4-h_5}}{\left( 1-x\left( 1-y\left( 1-z \right) \right)\right)^{h_1}}\,dxdydz\\ &=\frac{\Gamma(h_0+1)\prod_{j=2}^4\Gamma(h_j)\prod_{j=1}^4\Gamma(h_0+1-h_j-h_{j+1})}{\prod_{j=1}^5\Gamma(h_0+1-h_j)}\times\\ &\times\,_7F_6\left( \left. \begin{matrix}h_0,1+h_0/2,h_1,h_2,h_3,h_4,h_5\\h_0/2,1+h_0-h_1,1+h_0-h_2,1+h_0-h_3,1+h_0-h_4,1+h_0-h_5\end{matrix}\right|1 \right) \end{align} To evaluate \begin{align} \mathcal{I}&=\int_{0}^{1}\frac{x\mathbf{K}^2\left ( x \right )}{\sqrt{1-x^{2}}}\mathrm{d}x\\ &=\int_{0}^{1}\mathbf{K'}^2\left ( y \right )\mathrm{d}y \end{align} where $\mathbf{K'}\left ( x \right )=\mathbf{K}\left ( \sqrt{1-x^2} \right )$ (note that the arguments of the elliptic integrals are the elliptical parameter - not the elliptic modulus), we use the two integral representations \begin{align} \mathbf{K}\left ( s\right )&=\int_0^1\frac{dt}{\sqrt{(1-t^2)(1-s^2t^2)}}\\ &=\frac{1}{2}\int_0^1\frac{dx}{\sqrt{x(1-x)(1-s^2x)}} \tag{A} \end{align} and also, from the latter, \begin{align} \mathbf{K'}\left ( a\right )&=\frac{1}{2}\int_0^1\frac{du}{\sqrt{u(1-u)(1-(1-a^2)u)}}\\ &=\frac{1}{2}\int_{a^2}^1\frac{dv}{\sqrt{v(1-v)(v-a^2)}}\tag{B} \end{align} obtained using the change of variable $v=1-(1-a^2)u$. Then, \begin{align} \mathcal{I}&\overset{\text{(B)}}{=}\frac{1}{2}\int_0^1da\mathbf{K}(\sqrt{1-a^2})\int_{a^2}^1\frac{dy}{\sqrt{y(1-y)(y-a^2)}}\\ &=\frac{1}{2}\int_0^1dy \int_0^{\sqrt{y}}\frac{\mathbf{K}(\sqrt{1-a^2})}{\sqrt{y(1-y)(y-a^2)}}\,da\\ &=\frac{1}{4}\int_0^1dy \int_0^1\frac{\mathbf{K}(\sqrt{1-yz})}{\sqrt{y(1-y)(1-z)}}\,dz\\ &\overset{\text{(A)}}{=}\frac{1}{8}\int_0^1dy \int_0^1\frac{dz}{\sqrt{y(1-y)(1-z)}}\int_0^1\frac{dx}{\sqrt{x(1-x)(1-x(1-yz))}}\\ &\overset{z\to 1-z}{=}\frac{1}{8}\iiint_{\left[ 0,1\right]^3}\frac{(xyz)^{-1/2}\left[ (1-x)(1-y)(1-z)\right]^{-1/2}}{\sqrt{1-x\left( 1-y\left( 1-z \right) \right)}}dxdydz \end{align} With $h_j=\frac{1}{2}$ for $j=0,1,2,3,4,5$, using the integral representation for the hypergeometric series given above \begin{align} \mathcal{I}&= \frac{1}{8}\Gamma\left(\frac{3}{2}\right)\left[\Gamma\left(\frac{1}{2} \right) \right]^7 \,_7F_6\left( \left. \begin{matrix}1/2,1/2,1/2,1/2,1/2,1/2,5/4\\1/4,1,1,1,1,1\end{matrix}\right|1 \right) \end{align} which is the expected result.

Not a complete answer, but some (hopefully interesting) considerations, too long to fit in a comment.

I prefer to denote the complete elliptic integral of the first kind by regarding its argument as the elliptic modulus rather than the elliptic parameter. According to this convention (which is also Mathematica's one) the integral of interest is $$ \int_{0}^{1}\frac{x K^2(x^2)}{\sqrt{1-x^2}}\,dx =\frac{1}{2}\int_{0}^{1}\frac{K^2(x)}{\sqrt{1-x}}\,dx$$ where twice the RHS can be regarded as the $L^2$ squared norm of $\frac{K(x)}{\sqrt[4]{1-x}}$ over $(0,1)$.
The Fourier-Legendre series expansion of $K(x)$ is fairly simple and the Fourier-Legendre series expansion of $\frac{1}{\sqrt[4]{1-x}}$ can be computed from Rodrigues formula. To recover the mentioned result, it should be enough to compute a convolution, then apply Parseval's identity. Or write the original integral as $$ \int_{0}^{1}\int_{0}^{\pi/2}\int_{0}^{\pi/2}\frac{x}{\sqrt{(1-x^2)(1-x^2\cos^2\theta)(1-x^2\cos^2\phi)}}\,d\theta\,d\phi\,dx $$ or $$ \iiint_{(0,\pi/2)^3}\frac{\sin\eta}{\sqrt{(1-\sin^2\eta\cos^2\theta)(1-\sin^2\eta^2\cos^2\phi)}}\,d\theta\,d\phi\,d\eta $$ then apply Fubini's theorem / Beuker-Calabi-Kolk-type / reverse hyperspherical coordinates substitutions.