Expressing the integral $\int_{0}^{1}\frac{\mathrm{d}x}{\sqrt{\left(1-x^3\right)\left(1-a^6x^3\right)}}$ in terms of elliptic integrals
General Procedure:
I will give an outline of a possible procedure (that I learned from a calculus book some time ago) to solve integrals of the form
$$\mathcal{J}=\int\frac{{\rm d}x}{\sqrt{(x^2+ax+b)(x^2+cx+d)}}$$
First, let $x=t+\lambda$.
\begin{align}
\mathcal{J}=\int\frac{{\rm d}t}{\sqrt{(t^2+(2\lambda+a)t+\lambda^2+\lambda a+b)(t^2+(2\lambda+c)t+\lambda^2+\lambda c+d)}}
\end{align}
Next, let $\lambda=\dfrac{d-b}{a-c}$, then make the substitution $\displaystyle t=\mu v$ where $\displaystyle \mu=\sqrt{\lambda^2+\frac{ad-bc}{a-c}}$.
$$\mathcal{J}=\frac{1}{\mu}\int\frac{{\rm d}v}{\sqrt{\left(v^2+Av+1\right)\left(v^2+Cv+1\right)}}$$
where $\displaystyle A=\frac{2\lambda+a}{\mu}$ and $\displaystyle C=\frac{2\lambda+c}{\mu}$. Follow this up with the substitution $\displaystyle v=\frac{1+w}{1-w}$.
$$\mathcal{J}=\frac{2}{\mu}\int\frac{{\rm d}w}{\sqrt{((2-A)w^2+(2+A))((2-C)w^2+(2+C))}}$$
By applying the appropriate trigonometric substitution, it can be seen that $\mathcal{J}$ can be reduced to incomplete elliptic integrals.
Parameters
I will define the following parameters to save space and typing.
\begin{align}
\alpha_1=&\ \frac{1}{\sqrt{\left(a+1\right)^2\left(3a^2+2\sqrt{3}\sqrt{a^4+a^2+1}+3\right)}}\\
\alpha_2=&\ \frac{1}{\sqrt{\left(a-1\right)^2\left(3a^2+2\sqrt{3}\sqrt{a^4+a^2+1}+3\right)}}\\
\beta_1=&\ \frac{\sqrt{3}a^2-\sqrt{a^4+a^2+1}+\sqrt{3}a+\sqrt{3}}{\sqrt{3}a^2+\sqrt{a^4+a^2+1}+\sqrt{3}a+\sqrt{3}}\\
\beta_2=&\ \frac{\sqrt{3}a^2-\sqrt{a^4+a^2+1}-\sqrt{3}a+\sqrt{3}}{\sqrt{3}a^2+\sqrt{a^4+a^2+1}-\sqrt{3}a+\sqrt{3}}\\
\xi_1=&\ \frac{\sqrt{a^4+a^2+1}-\sqrt{3}a}{\sqrt{a^4+a^2+1}+\sqrt{3}a}\\
\xi_2=&\ \frac{\sqrt{a^4+a^2+1}+\sqrt{3}a}{\sqrt{a^4+a^2+1}-\sqrt{3}a}
\end{align}
Evaluation of $\ \mathcal{I}(a)$
It was derived by the OP that
\begin{align}
\mathcal{I}(a)
=&\underbrace{\frac{1}{4a}\int^\infty_\frac{a^2+1}{2a}\frac{{\rm d}x}{\sqrt{\left(x^2-\frac{(a-1)^2}{2a}x-\frac{a^2+1}{2a}\right)\left(x^2+\frac{a^2+1}{2a}x+\frac{a^4-a^2+1}{4a^2}\right)}}}_{\mathcal{I}_1}\\
&+\underbrace{\frac{1}{4a}\int^\infty_\frac{a^2+1}{2a}\frac{{\rm d}x}{\sqrt{\left(x^2-\frac{(a+1)^2}{2a}x+\frac{a^2+1}{2a}\right)\left(x^2+\frac{a^2+1}{2a}x+\frac{a^4-a^2+1}{4a^2}\right)}}}_{\mathcal{I}_2}
\end{align}
Applying the aforementioned procedure and using Wolfram Alpha to do all the algebra, we get
\begin{align}
\mathcal{I}_1
=&\frac{1}{4a}\int^\infty_\frac{a^2+1}{2a}\frac{{\rm d}x}{\sqrt{\left(x^2-\frac{(a-1)^2}{2a}x-\frac{a^2+1}{2a}\right)\left(x^2+\frac{a^2+1}{2a}x+\frac{a^4-a^2+1}{4a^2}\right)}}\\
=&\frac{1}{\sqrt{3}\sqrt{a^4+a^2+1}}\int^\infty_\frac{\sqrt{3}(a^2+a+1)}{\sqrt{a^4+a^2+1}}\frac{{\rm d}x}{\sqrt{\left(x^2-\frac{2(2a^2+a+2)}{\sqrt{3}\sqrt{a^4+a^2+1}}x+1\right)\left(x^2-\frac{2\sqrt{3}\ a}{\sqrt{a^4+a^2+1}}+1\right)}}\\
=&\alpha_1\int^1_{\beta_1}\frac{{\rm d}x}{\sqrt{\left(x^2-\beta_1^2\right)\left(x^2+\xi_1\right)}}
=\alpha_1\int^{\arccos{\beta_1}}_{0}\frac{\sec{\varphi}\tan{\varphi}}{\sqrt{(\sec^2{\varphi}-1)(\beta_1^2\sec^2{\varphi}+\xi_1)}}{\rm d}\varphi\\
=&\alpha_1\int^{\arccos{\beta_1}}_{0}\frac{{\rm d}\varphi}{\sqrt{\beta_1^2+\xi_1\cos^2{\varphi}}}
=\alpha_1\int^{\arccos{\beta_1}}_{0}\frac{{\rm d}\varphi}{\sqrt{(\beta_1^2+\xi_1)-\xi_1\sin^2{\varphi}}}\\
=&\frac{\alpha_1}{\sqrt{\beta_1^2+\xi_1}}\int^{\arccos{\beta_1}}_{0}\frac{{\rm d}\varphi}{\sqrt{1-\frac{\xi_1}{\beta_1^2+\xi_1}\sin^2{\varphi}}}
=\color{red}{\frac{\alpha_1}{\sqrt{\beta_1^2+\xi_1}}\mathbf{F}\left(\arccos{\beta_1},\sqrt{\frac{\xi_1}{\beta_1^2+\xi_1}}\right)}
\end{align}
for $\mathcal{I}_1$, and we obtain (I omit the steps of trigonometric manipulations since they are the same)
\begin{align}
\mathcal{I}_2
=&\frac{1}{4a}\int^\infty_\frac{a^2+1}{2a}\frac{{\rm d}x}{\sqrt{\left(x^2-\frac{(a+1)^2}{2a}x+\frac{a^2+1}{2a}\right)\left(x^2+\frac{a^2+1}{2a}x+\frac{a^4-a^2+1}{4a^2}\right)}}\\
=&\frac{1}{\sqrt{3}\sqrt{a^4+a^2+1}}\int^\infty_{\frac{\sqrt{3}(a^2-a+1)}{\sqrt{a^4+a^2+1}}}\frac{{\rm d}x}{\sqrt{\left(x^2-\frac{2(2a^2-a+2)}{\sqrt{3}\sqrt{a^4+a^2+1}}x+1\right)\left(x^2+\frac{2\sqrt{3}\ a}{\sqrt{a^4+a^2+1}}x+1\right)}}\\
=&\alpha_2\int^1_{\beta_2}\frac{{\rm d}x}{\sqrt{\left(x^2-\beta_2^2\right)\left(x^2+\xi_2\right)}}=\color{blue}{\frac{\alpha_2}{\sqrt{\beta_2^2+\xi_2}}\mathbf{F}\left(\arccos{\beta_2},\sqrt{\frac{\xi_2}{\beta_2^2+\xi_2}}\right)}
\end{align}
Therefore
$$\mathcal{I}(a)=\frac{\alpha_1}{\sqrt{\beta_1^2+\xi_1}}\mathbf{F}\left(\arccos{\beta_1},\sqrt{\frac{\xi_1}{\beta_1^2+\xi_1}}\right)+\frac{\alpha_2}{\sqrt{\beta_2^2+\xi_2}}\mathbf{F}\left(\arccos{\beta_2},\sqrt{\frac{\xi_2}{\beta_2^2+\xi_2}}\right)$$
I have substituted in some random values of $a$ into $\mathcal{I}_1$ and it seems to match numerically. However, I have not verified $\mathcal{I}_2$ yet so there might be a mistake.