Lemma 1: For $\lambda \in \mathbb{R}^{+}$ $$\sum\limits_{n=1}^{\infty} \dfrac{\lambda^n}{n!}\mathcal{H}_{n} = e^{\lambda}\left(\ln \lambda + \gamma - \operatorname{Ei}(-\lambda)\right) \tag{1} \label{lemma1}$$
where $\mathcal{H}_{n}$ — $n$-th harmonic number and $\operatorname{Ei}(\cdot)$ — exponential integral.
Start from two series representations of lower incomplete gamma function $\gamma(\beta, \lambda)$:
\begin{align}
\gamma(\beta, \lambda) &= e^{-\lambda}\sum\limits_{n=0}^{\infty} \dfrac{\lambda^{n+\beta}}{\beta\left(\beta+1\right)\ldots\left(\beta+n\right)} \\
\gamma(\beta, \lambda) &= \sum\limits_{n=0}^{\infty} (-1)^{n}\dfrac{\lambda^{n+\beta}}{n!\left(\beta+n\right)}
\end{align}
Now take the derivative with respect to $\beta$ at $1$. Since
\begin{align*}
\dfrac{\mathrm{d}}{\mathrm{d}\beta}\left(\dfrac{1}{\beta\left(\beta+1\right)\ldots\left(\beta+n\right)}\right)_{\beta=1} &= -\left.\dfrac{1}{\beta\left(\beta+1\right)\ldots\left(\beta+n\right)}\left(\dfrac{1}{\beta}+\dfrac{1}{\beta+1}+\ldots+\dfrac{1}{\beta+n}\right)\right\vert_{\ \beta=1} \\
&= -\dfrac{\mathcal{H}_{n+1}}{(n+1)!}
\end{align*}
we have that
\begin{align*}
-e^{-\lambda}\sum\limits_{n=1}^{\infty} \dfrac{\lambda^{n}}{n!}\mathcal{H}_{n}+\ln \lambda\ e^{-\lambda}\sum\limits_{n=1}^{\infty} \dfrac{\lambda^{n}}{n!} &= -\ln \lambda\sum\limits_{n=1}^{\infty} (-1)^{n}\dfrac{\lambda^{n}}{n!} + \sum\limits_{n=1}^{\infty} (-1)^{n}\dfrac{\lambda^{n}}{n!\ n} \\
-e^{-\lambda}\sum\limits_{n=1}^{\infty} \dfrac{\lambda^{n}}{n!}\mathcal{H}_{n}+\ln \lambda \left(1-e^{-\lambda}\right) &= \ln \lambda \left(1-e^{-\lambda}\right) + \sum\limits_{n=1}^{\infty} (-1)^{n}\dfrac{\lambda^{n}}{n!\ n} \\
\sum\limits_{n=1}^{\infty} \dfrac{\lambda^{n}}{n!}\mathcal{H}_{n} &= -e^{\lambda}\sum\limits_{n=1}^{\infty} (-1)^{n}\dfrac{\lambda^{n}}{n!\ n} \\
&= e^{\lambda}\left(\ln \lambda + \gamma - \operatorname{Ei}(-\lambda)\right)
\end{align*}
Let
$$
\mathfrak{I}(\lambda, \alpha, \beta) = \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}\,x^{\alpha - 1}}{\left(1+x\right)^{\alpha+\beta}}\,\mathrm{d}x
$$
Lemma 2: Let $\alpha, \lambda \in \mathbb{R}^{+}$ and $\alpha + \beta > 0$. Then
$$
\mathfrak{I}(\lambda, \alpha, \beta) = \dfrac{\Gamma(\alpha)}{\Gamma(\alpha+\beta)}\lambda^{\beta}\mathfrak{I}(\lambda, \alpha + \beta, -\beta) \tag{2} \label{lemma2}
$$
Using Laplace transform properties we have
\begin{align*}
\mathfrak{I}(\lambda, \alpha, \beta) &= \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}x^{\alpha - 1}}{\left(1+x\right)^{\alpha+\beta}}\,\mathrm{d}x \\
&= \int\limits_{0}^{\infty} e^{-\lambda t}t^{\alpha-1}\,\mathcal{L}\left\{\dfrac{e^{-x}x^{\alpha+\beta-1}}{\Gamma(\alpha+\beta)}\right\}(t)\,\mathrm{d}t \\
&= \int\limits_{0}^{\infty} \dfrac{\Gamma(\alpha)}{\left(\lambda+x\right)^{\alpha}}\dfrac{e^{-x}x^{\alpha+\beta-1}}{\Gamma(\alpha+\beta)} \,\mathrm{d}x \\
&= \dfrac{\Gamma(\alpha)}{\Gamma(\alpha+\beta)}\lambda^{\beta}\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}x^{\alpha+\beta-1}}{\left(1+x\right)^{\alpha}}\,\mathrm{d}x \\
&= \dfrac{\Gamma(\alpha)}{\Gamma(\alpha+\beta)}\lambda^{\beta}\mathfrak{I}(\lambda, \alpha + \beta, -\beta)
\end{align*}
Claim 1: For $\lambda > 0$ and $\beta \in \mathbb{R}$
$$
\mathfrak{I}(\lambda, 1, \beta) = e^{\lambda}\lambda^{\beta}\Gamma(-\beta, \lambda) \tag{3} \label{claim1}
$$
where $\Gamma(\cdot, \cdot)$ — upper incomplete gamma function.
\begin{align*}
\mathfrak{I}(\lambda, 1, \beta) &= \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{\left(1+x\right)^{1+\beta}}\,\mathrm{d}x \\
&= e^{\lambda}\int\limits_{1}^{\infty} \dfrac{e^{-\lambda x}}{x^{1+\beta}}\,\mathrm{d}x \\
&= e^{\lambda}\lambda^{\beta}\int\limits_{\lambda}^{\infty} \dfrac{e^{-x}}{x^{1+\beta}}\,\mathrm{d}x \\
&= e^{\lambda}\lambda^{\beta}\Gamma(-\beta, \lambda)
\end{align*}
After applying \eqref{lemma2} to \eqref{claim1} we have that
$$
\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}x^{\beta}}{1+x}\,\mathrm{d}x = e^{\lambda}\Gamma(1+\beta)\Gamma(-\beta, \lambda) \tag{4} \label{consec1}
$$
Taking the derivative with respect to the variable $\beta$ at point $(\lambda, \beta) = (1, 0)$ leads us to
\begin{align*}
\int\limits_{0}^{\infty} \dfrac{e^{-x}\ln x}{1+x}\,\mathrm{d}x &= e\dfrac{\partial}{\partial \beta}\Big(\Gamma(1+\beta)\Gamma(-\beta, 1)\Big)\Bigg\vert_{\beta=0} \\
&= -e\gamma \int\limits_{1}^{\infty} \dfrac{e^{-x}}{x}\,\mathrm{d}x -e \int\limits_{1}^{\infty} \dfrac{e^{-x}\ln x}{x}\,\mathrm{d}x \\
&= e\gamma\operatorname{Ei}(-1)-\dfrac{1}{2}e\int\limits_{1}^{\infty} e^{-x}\ln^2 x\,\mathrm{d}x \\
&= e\gamma\operatorname{Ei}(-1)-\dfrac{1}{2}e\dfrac{\mathrm{d}^2}{\mathrm{d}x^2}\Gamma(1)+\dfrac{1}{2}e\int\limits_{0}^{1} e^{-x}\ln^2 x\,\mathrm{d}x \\
\end{align*}
$$
\int\limits_{0}^{\infty} \dfrac{e^{-x}\ln x}{1+x}\,\mathrm{d}x = e\gamma\operatorname{Ei}(-1)-\dfrac{1}{2}e\left(\gamma^2+\dfrac{1}{6}\pi^2\right)-e\sum\limits_{n=1}^{\infty} \dfrac{(-1)^n}{n!\,n^2} \tag{5} \label{first}
$$
Let
$$
f(\lambda) = \int\limits_{0}^{\infty} e^{-x}\ln^2 \left(\lambda x + x^2\right)\,\mathrm{d}x
$$
Then
\begin{align*}
f'(\lambda) &= 2\int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(\lambda x + x^2\right)}{\lambda + x}\,\mathrm{d}x \\
&= 2\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{1+x}\Big(\ln \left(1+x\right)+\ln x + 2 \ln \lambda\Big)\,\mathrm{d}x \\
&= 4\ln \lambda \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{1+x}\,\mathrm{d}x + 2\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{1+x}\ln \left(1+x\right)\,\mathrm{d}x + 2\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{1+x}\ln x\,\mathrm{d}x \\
&= -4\ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda) - 2\dfrac{\partial}{\partial \beta}\left.\left(\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{\left(1+x\right)^{1+\beta}}\,\mathrm{d}x\right)\right\vert_{\beta=0} + 2\dfrac{\partial}{\partial \beta}\left.\left(\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}x^{\beta}}{1+x}\,\mathrm{d}x\right)\right\vert_{\beta=0} \\
&= -4\ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda) -2e^{\lambda}\dfrac{\partial}{\partial \beta}\Big(\lambda^{\beta}\Gamma(-\beta,\lambda)\Big)\Bigg\vert_{\beta=0} + 2e^{\lambda}\dfrac{\partial}{\partial \beta}\Big(\Gamma(1+\beta)\Gamma(-\beta,\lambda)\Big)\Bigg\vert_{\beta=0} \\
&= -2\ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda)+2\gamma e^{\lambda}\operatorname{Ei}(-\lambda)
\end{align*}
Integrating over $[0,1]$ leads us to
\begin{align*}
f(1) &= f(0) -2\int\limits_{0}^{1} \ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda)\,\mathrm{d}\lambda + 2\gamma\Big(e^{\lambda}\operatorname{Ei}(-\lambda)-\ln \lambda\Big)\Bigg\vert_{0}^{1} \\
&= 2\gamma e\operatorname{Ei}(-1)+2\gamma^2+\dfrac{2}{3}\pi^2-2\int\limits_{0}^{1} \ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda)\,\mathrm{d}\lambda
\end{align*}
For integral in last formula we use \eqref{lemma1}:
\begin{align*}
\int\limits_{0}^{1} \ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda)\,\mathrm{d}\lambda &= \int\limits_{0}^{1} \ln^2 \lambda e^{\lambda}\,\mathrm{d}\lambda +\gamma \int\limits_{0}^{1} \ln \lambda e^{\lambda}\,\mathrm{d}\lambda + \sum\limits_{n=2}^{\infty} \dfrac{\mathcal{H}_{n-1}}{n!\ n} \\
&= \sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2} + \gamma^2-\gamma\operatorname{Ei}(1) + \sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n}
\end{align*}
Combining all together we have
$$
\int\limits_{0}^{\infty} e^{-x}\ln^2 \left(x+x^2\right)\,\mathrm{d}x = 2\gamma\Big(e\operatorname{Ei}(-1)+\operatorname{Ei}(1)\Big)+\dfrac{2}{3}\pi^2-2\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}-2\sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n} \tag{6} \label{second}
$$
From another side
\begin{align*}
\int\limits_{0}^{\infty} e^{-x}\ln^2 \left(x+x^2\right)\,\mathrm{d}x &= \int\limits_{0}^{\infty} e^{-x}\Big(\ln^2 \left(1+x\right) + \ln^2 x + 2\ln \left(1+x\right)\ln x\Big)\,\mathrm{d}x \\
&= \gamma^2+\dfrac{1}{6}\pi^2+e\int\limits_{1}^{\infty}e^{-x}\ln^2 x\,\mathrm{d}x + 2\int\limits_{0}^{\infty} e^{-x}\ln \left(1+x\right)\ln x\,\mathrm{d}x \\
&= \left(1+e\right)\left(\gamma^2+\dfrac{1}{6}\pi^2\right)+2e\sum\limits_{n=1}^{\infty} \dfrac{(-1)^n}{n!\ n^2} + 2\int\limits_{0}^{\infty} e^{-x}\ln \left(1+x\right)\ln x\,\mathrm{d}x
\end{align*}
Express last integral
\begin{align*}
\int\limits_{0}^{\infty} e^{-x}\ln \left(1+x\right)\ln x\,\mathrm{d}x &= \gamma\Big(e\operatorname{Ei}(-1)+\operatorname{Ei}(1)\Big) -\dfrac{1}{2}\gamma^2+\dfrac{1}{4}\pi^2-\dfrac{1}{2}e\left(\gamma^2+\dfrac{1}{6}\pi^2\right) \\
&-e\sum\limits_{n=1}^{\infty} \dfrac{(-1)^n}{n!\ n^2}-\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}-\sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n}
\end{align*}
Now use integration by parts and already obtained result \eqref{first}:
\begin{align*}
\int\limits_{0}^{\infty} e^{-x}\ln \left(1+x\right)\ln x\,\mathrm{d}x &= \int\limits_{0}^{\infty} e^{-x}\dfrac{\ln x}{1+x}\,\mathrm{d}x + \int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(1+x\right)}{x}\,\mathrm{d}x \\
&= e\gamma\operatorname{Ei}(-1)-\dfrac{1}{2}e\left(\gamma^2+\dfrac{1}{12}\pi^2\right)-e\sum\limits_{n=1}^{\infty} \dfrac{(-1)^n}{n!\ n^2}+ \int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(1+x\right)}{x}\,\mathrm{d}x
\end{align*}
So
$$
\int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(1+x\right)}{x}\,\mathrm{d}x = -\dfrac{1}{2}\gamma^2+\dfrac{1}{4}\pi^2+\gamma\operatorname{Ei}(1)-\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}-\sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n} \tag{7} \label{final}
$$
Back to original problem
With integration by parts and substitution original integral can be converted to
\begin{align*}
\int\limits_{0}^{e} \operatorname{Li}_2\left(\ln x\right)\,\mathrm{d}x &= \dfrac{1}{6}e\pi^2-\int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(1+x\right)}{x}\,\mathrm{d}x+\int\limits_{0}^{1} e^{x}\dfrac{\ln \left(1-x\right)}{x}\,\mathrm{d}x \\
&= \dfrac{1}{6}e\pi^2 + \dfrac{1}{2}\gamma^2-\dfrac{1}{4}\pi^2-\gamma\operatorname{Ei}(1)+\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}+\sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n} \\
&+\int\limits_{0}^{1} \dfrac{\ln\left(1-x\right)}{x}\,\mathrm{d}x+\sum\limits_{n=1}^{\infty}\dfrac{1}{n!}\int\limits_{0}^{1} x^{n-1}\ln\left(1-x\right)\,\mathrm{d}x \\
&= \dfrac{1}{6}e\pi^2 + \dfrac{1}{2}\gamma^2-\dfrac{5}{12}\pi^2-\gamma\operatorname{Ei}(1)+\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}
\end{align*}
