Partial answer. Regarding $I$ as a function of $u$, we note that
\begin{align*}
I(0) &= \int_{0}^{\infty} x^2 \log(1-e^{-x})\, \mathrm{d}x = -\frac{\pi^4}{45}, \\
I'(0) &= 0, \\
I''(0) &= \int_{0}^{\infty} \frac{x}{e^x - 1} \, \mathrm{d}x = \frac{\pi^2}{6}.
\end{align*}
To get the higher terms, write
\begin{align*}
I(u)
&= I(0^+) + \frac{I''(0)}{2}u^2 + \int_{0}^{\infty} \bigg( \underbrace{ x^2 \log\left(1 - e^{-\sqrt{x^2+u^2}}\right) - x^2 \log\left(1 - e^{-x}\right) - \frac{u^2 x}{2(e^x - 1)} }_{=\text{(*)}} \bigg) \, \mathrm{d}x.
\end{align*}
Here, the integrand $\text{(*)}$ can be recast as
\begin{align*}
\text{(*)}
&= \int_{0}^{u} \left[ \frac{\partial}{\partial s} x^2 \log\left(1 - e^{-\sqrt{x^2+u^2}}\right) - \frac{sx}{(e^x - 1)} \right]\, \mathrm{d}s \\
&= \int_{0}^{u} \bigg[ \frac{sx^2}{\sqrt{x^2+s^2}(e^{\sqrt{x^2+s^2}}-1)} - \frac{sx}{(e^x - 1)} \bigg]\, \mathrm{d}s \\
&= \int_{0}^{u} \int_{0}^{s} s x^2 \frac{\partial}{\partial t} \bigg( \frac{1}{\sqrt{x^2+t^2}(e^{\sqrt{x^2+t^2}}-1)} \bigg) \, \mathrm{d}t \mathrm{d}s \\
&= -\int_{0}^{u} \int_{0}^{s} \frac{s x^2 t}{(t^2 + x^2)^2} f\left(\sqrt{x^2+t^2}\right) \, \mathrm{d}t\mathrm{d}s,
\end{align*}
where $f(a) = \frac{a(e^a(1+a) - 1)}{(e^a - 1)^2}$. For the future use, we remark that $a = 0$ is a removable singularity of $f$ with $f(0) = 2$ and $2 \geq f(a) \geq f(b) \geq 0$ for all $0 \leq a \leq b$. Interchanging the order of integration,
\begin{align*}
\text{(*)}
&= -\int_{0}^{u} \frac{x^2 t (u^2-t^2)}{2(t^2 + x^2)^2} f\left(\sqrt{x^2+t^2}\right) \, \mathrm{d}t.
\end{align*}
Plugging this back and applying the substitution $(x, t) \mapsto (utx, ut)$,
\begin{align*}
I(u)
&= I(0^+) + \frac{I''(0)}{2}u^2 - u^3 \int_{0}^{\infty} \int_{0}^{1} \frac{x^2 (1-t^2)}{2(1 + x^2)^2} f\left(u t \sqrt{x^2+1}\right) \, \mathrm{d}t\mathrm{d}x.
\end{align*}
As $u \to 0^+$, monotone convergence tells that
\begin{align*}
&\int_{0}^{\infty} \int_{0}^{1} \frac{x^2 (1-t^2)}{2(1 + x^2)^2} f\left(u t \sqrt{x^2+1}\right) \, \mathrm{d}t\mathrm{d}x \\
&\xrightarrow[u \to 0^+]{} \int_{0}^{\infty} \int_{0}^{1} \frac{x^2 (1-t^2)}{(1 + x^2)^2} \, \mathrm{d}t\mathrm{d}x
= \frac{\pi}{6}.
\end{align*}
So it follows that
$$ I(u) = -\frac{\pi^4}{45} + \frac{\pi^2}{12}u^2 - \frac{\pi}{6}u^3 + o(u^3). $$
Numerical analysis. If $0 < x < \frac{1}{2}$, then $\left| \log (1 - x) \right| = -\log(1-x) \leq \frac{x}{1-x} \leq 2x $. Using this, we note that, for $R > 1$,
$$ \left|\int_{R}^{\infty} x^2 \log\left( 1 - e^{-\sqrt{x^2 + u^2}} \right) \, \mathrm{d}x\right|
\leq \int_{R}^{\infty} 2x^2 e^{-x} \, \mathrm{d}x
= 2(R^2 + 2R + 2) e^{-R}. $$
For instance, this bound is $< 10^{-30}$ for $R \geq 80$, and so, we may compute
$$ \int_{0}^{\infty} x^2 \log\left( 1 - e^{-\sqrt{x^2 + u^2}} \right) \, \mathrm{d}x = \int_{0}^{80} x^2 \log\left( 1 - e^{-\sqrt{x^2 + u^2}} \right) \, \mathrm{d}x \pm 10^{-30}$$
uniformly in $u > 0$. Since the integrand is continuous on $[0, 80]$, we may use any CAS to easily estimate the truncated integral with a precision of $10^{-30}$. Then for the quantities
\begin{align*}
\text{(A)} &= \text{numerical integration of } \int_{0}^{80} x^2 \log\left( 1 - e^{-\sqrt{x^2 + u^2}} \right) \, \mathrm{d}x, \\
\text{(B)} &= \text{OP's asymtotic form} \\
&= -\frac{\pi^4}{45} + \frac{\pi^2}{12}u^2 - \frac{\pi}{6}u^3 + \frac{\frac{3}{2} + 2\log(4\pi) - 2\gamma - 2\log u}{32} u^4, \\
\text{(C)} &= \text{@Yuri Negometyanov's asymptotic form} \\
&= -\frac{u^2}{2}\operatorname{Li}_2(e^{-u}) -2u \operatorname{Li}_3(e^{-u}) - 2\operatorname{Li}_4(e^{-u}),
\end{align*}
we obtain the following comparisons of numerical values:

Coinciding leading decimals are highlighted as red.