integration - Closed form for $int_0^emathrm{Li}_2(ln{x}),dx$?

Inspired by this question and this answer, I decided to investigate the family of integrals
$$I(k)=\int_0^e\mathrm{Li}_k(\ln{x})\,dx,\tag{1}$$
where $\mathrm{Li}_k(z)$ represents the polylogarithm of order $k$ and argument $z$. $I(1)$ evaluates to $e\gamma$, but $I(2)$ has resisted my efforts (which can be seen here).

Neither ISC nor WolframAlpha could provide a closed form for its numerical value--however, I've conjectured a possible analytic form. $$\eqalign{\int_0^e\mathrm{Li}_2(\ln{x})\,dx&\stackrel?=\,_3F_3(1,1,1;2,2,2;1)+\frac{\pi^2(2e-5)}{12}+\frac{\gamma^2}{2}-\gamma\,\mathrm{Ei}(1)\\&=0.578255559804073275225659054377625577...\tag{2}}$$

Brevan Ellefsen has computed that my conjecture is accurate to at least 150 digits. Brevan also gave the alternate form
$$\frac{\pi^2e}{6}+\gamma G_{1,2}^{2,0}\left(-1\left|\begin{array}{c}1\\0,0\\\end{array}\right.\right)+G_{2,3}^{3,0}\left(-1\left|\begin{array}{c}1,1\\0,0,0\\\end{array}\right.\right).\tag{2.1}$$

Is there a closed form for $I(2)$ that doesn't involve Meijer G or hypergeometrics? The simplicity of the following two equations seems to suggest that there might be. $(3.1)$ follows directly from $(3)$, which I've proven here.

$$\begin{align} \sum_{k=1}^\infty I(k)&=e\tag{3}\\\sum_{k=2}^\infty I(k)&=e(1-\gamma)\tag{3.1}\end{align}$$

---

PROGRESS UPDATE: Using this equation, I've turned $_3F_3(1,1,1;2,2,2;1)$ into
$$\lim_{c\to 1}\left(\frac{\mathrm{Ei}(1)-\gamma}{c-1}+\frac{1-e}{(c-1)^2}+\frac{(-1)^{-c}\,\Gamma(c-1)}{c-1}+\frac{(-1)^{1-c}\,\Gamma(c,-1)}{(c-1)^2}\right),\tag{4}\label{4}$$ but I don't know how to proceed from there. EDIT: This limit leads nowhere. See below.

SECOND PROGRESS UPDATE: After some studying of the properties of the Meijer G function, I've finally cracked the limit; however, the result is an underwhelming $_3F_3(1,1,1;2,2,2;1)$. Before I evaluate the limit, I'd first like to state the following intermediate result:

Lemma $(4.1)$: For $z\in\mathbb{C}$,
$$G_{2,3}^{3,0}\left(z\left|\begin{array}{c}1,1\\0,0,0\\\end{array}\right.\right)=\gamma\ln{z}+\frac12\ln^2(z)-z\,_3F_3(1,1,1;2,2,2;-z)+\frac{\gamma^2}2+\frac{\pi^2}{12}.\tag{4.1}\label{4.1}$$

My proof for this can be found here. Now I return to the limit $\eqref{4}$.
Consider the following: $\frac{1}{c-1}=\frac{c-1}{(c-1)^2}$, $(c-1)\Gamma(c-1)=\Gamma(c)$, and $(-1)^{1-c}=-(-1)^{-c}$. Based on these algebraic identities, the limit can be written as
$$\lim_{c\to 1}\frac{(c-1)(\mathrm{Ei}(1)-\gamma)+1-e+(-1)^{-c}(\Gamma(c)-\Gamma(c,-1))}{(c-1)^2}.\tag{4.2}$$
In this form, the limit is $\frac{0}0$. Using l'Hospital twice, we obtain

$$\lim_{c\to 1}{(-1)^{-c}\left(-G_{3,4}^{4,0}\left(-1\left|\begin{array}{c}1,1,1\\0,0,0,c\\\end{array}\right.\right)+\Gamma{(c)}\left(\frac{\psi_0{(c)}^2}2-i\pi\psi_0(c)+\frac{\psi_1(c)}2-\frac{\pi^2}2\right)\right)}$$ $$\begin{align}&=G_{3,4}^{4,0}\left(-1\left|\begin{array}{c}1,1,1\\0,0,0,1\\\end{array}\right.\right)-\frac{\psi_0(1)^2}2+i\pi\psi_0(1)-\frac{\psi_1(1)}2+\frac{\pi^2}2\\&=G_{2,3}^{3,0}\left(-1\left|\begin{array}{c}1,1\\0,0,0\\\end{array} \right.\right)-\frac{\gamma^2}{2}-i\gamma\pi+\frac{5\pi^2}{12}.\tag{4.2a}\end{align}$$
Using Lemma $\eqref{4.1}$, we know that
$$G_{2,3}^{3,0}\left(-1\left|\begin{array}{c}1,1\\0,0,0\\\end{array}\right.\right)={}_3F_3(1,1,1;2,2,2;1)+\frac{\gamma^2}2+i\gamma\pi-\frac{5\pi^2}{12},\tag{4.3}$$ which can be rewritten as
$$G_{2,3}^{3,0}\left(-1\left|\begin{array}{c}1,1\\0,0,0\\\end{array}\right.\right)-\frac{\gamma^2}{2}-i\gamma\pi+\frac{5\pi^2}{12}={}_3F_3(1,1,1;2,2,2;1).\tag{4.3a}$$

Thus,
$$\eqalign{&\lim_{c\to 1}\left(\frac{\mathrm{Ei}(1)-\gamma}{c-1}+\frac{1-e}{(c-1)^2}+\frac{(-1)^{-c}\,\Gamma(c-1)}{c-1}+\frac{(-1)^{1-c}\,\Gamma(c,-1)}{(c-1)^2}\right)\\=&{}_3F_3(1,1,1;2,2,2;1).\tag{4.4}}$$