Here is an algorithmic way to get the asymptotic series (in powers of $1/a$ as $a\to\infty$).
The integral (after $x=t^{1/a}$) equals $f(1/a)$, where
$$
f(\lambda)=\int_0^1\frac{\lambda t^{\lambda-1}}{(1+t^\lambda+(1-t)^\lambda)^2}\,dt=\frac16+\lambda\int_0^1 F(\log t, \log(1-t),\lambda)\frac{dt}t,
$$
and
$$
F(a,b,\lambda)=\frac{e^{a\lambda}}{(1+e^{a\lambda}+e^{b\lambda})^2}-\frac{e^{a\lambda}}{(2+e^{a\lambda})^2}
$$
has a power series (that can be computed) of the form
$$
F(a,b,\lambda)=\sum_{n=1}^\infty F_n(a,b)\lambda^n,\qquad F_n(a,b)=\sum_{k=1}^n F_{n,k}a^{n-k}b^k.
$$
This gives the asymptotics $\color{blue}{f(\lambda)\asymp 1/6+\sum_{n=1}^\infty f_n\lambda^{n+1}}$ as $\lambda\to 0^+$, where
$$
f_n=\sum_{k=1}^n F_{n,k}G_{n,k},\quad G_{n,k}=L_{n-k,k},\quad L_{m,n}=\int_0^1\log^m(t)\log^n(1-t)\frac{dt}t.
$$
An algorithmic way to express $L_{m,n}$ in terms of $\zeta$-function is to consider
\begin{align*}
1+\alpha\sum_{m=0}^\infty\sum_{n=1}^\infty L_{m,n}\frac{\alpha^m}{m!}\frac{\beta^n}{n!}
&=1+\alpha\int_0^1 e^{\alpha\log t}(e^{\beta\log(1-t)}-1)\frac{dt}t
\\&=\alpha\int_0^1 t^{\alpha-1}(1-t)^\beta\,dt
=\frac{\Gamma(1+\alpha)\Gamma(1+\beta)}{\Gamma(1+\alpha+\beta)}
\end{align*}
and use known series $\log\Gamma(1+z)=-\gamma z+\sum_{n=2}^\infty\zeta(n)(-z)^n/n$; we obtain
$$
1+\alpha\sum_{m=0}^\infty\sum_{n=1}^\infty L_{m,n}\frac{\alpha^m}{m!}\frac{\beta^n}{n!}=\exp\left(\sum_{n=2}^\infty\frac{(-1)^n}n\zeta(n)\big[\alpha^n+\beta^n-(\alpha+\beta)^n\big]\right).$$
The exponentiation has to be performed algorithmically.
I've implemented the above in PARI/GP as follows:
experiment(N) = {
my (Z = vector(N, n, if (n % 2, eval(Str("Z", n)),
(2 * 'pi)^n / 2 / n! * abs(bernfrac(n)))));
Z = sum(n = 2, N, Z[n] * (-x)^n * (1 + y^n - (1 + y)^n) / n);
Z = apply(Vecrev, Vec((exp(Z + O(x^(N + 1))) - 1) / (x * y)));
my (F = exp(x) + O(x^N));
F = F / (1 + F + exp(x * y))^2 - F / (2 + F)^2;
F = apply(Vecrev, Vec(F / (x * y)));
return (vector(N - 1, n, sum(k = 1, #F[n],
k! * (n - k)! * F[n][k] * Z[n][k])))
};
so that running experiment(8) produces
\begin{align*}
f_1&=\frac{\pi^2}{81},
\\ f_2 &= 0,
\\ f_3 &= -\frac{\pi^4}{972},
\\ f_4 &= \frac{2\pi^2}{729}\zeta(3),
\\ f_5 &= \frac{7\pi^6}{29160},
\\ f_6 &= -\frac{5\pi^4}{4374}\zeta(3)-\frac{10\pi^2}{729}\zeta(5),
\\ f_7 &= -\frac{38047\pi^8}{396809280}-\frac{14\pi^2}{6561}\zeta^2(3),
\end{align*}
and so on. The values of $f_1$ and $f_3$ (and $f_2$ of course) are the same as reported by @Svyatoslav.
But, as we see, the asymptotics still has non-zero terms with odd powers of $\lambda$ (thus of $a$).
I have no idea whether a closed form of this asymptotic expansion exists.
