Contour integration with infinite poles

Let's assume, for the moment, that $k\leq 0$. Then we can rewrite our integral in the following way:

$$\int_{-\infty}^\infty \frac{x^2}{e^k + e^{x^2}}dx = e^{-k}\int_{-\infty}^\infty \frac{x^2e^ke^{-x^2}}{e^ke^{-x^2}+1}dx = e^{-k}\sum_{n=1}^\infty (-1)^{n+1}e^{nk}\int_{-\infty}^\infty x^2e^{-nx^2}dx$$

Then, by Feynman's trick

$$= e^{-k}\sum_{n=1}^\infty (-1)^{n+1} e^{nk}\Biggr(-\frac{d}{dn}\int_{-\infty}^\infty e^{-nx^2}dx\Biggr) = e^{-k}\sum_{n=1}^\infty (-1)^{n+1} e^{nk} \left(-\frac{d}{dn}\sqrt{\frac{\pi}{n}}\right)$$

$$ = -\frac{\sqrt{\pi}}{2}e^{-k} \sum_{n=1}^\infty \frac{(-e^k)^n}{n^{\frac{3}{2}}} \equiv -\frac{\sqrt{\pi}}{2}e^{-k} \text{Li}_{\frac{3}{2}}(-e^k)$$

While these specific steps don't apply when $k>0$, we can use analytic continuation to extend the polylogarithm (and taking the antilimit).

Using a CAS, I did not find any antiderivative. However, for the integral $$\int_{-\infty}^{\infty} \frac{x^2}{K+e^{x^2}}\,dx=-\frac{\sqrt{\pi } }{2 K} \text{Li}_{\frac{3}{2}}(-K)$$ where appears the polylogarithm function.

$$\int_{-\infty}^{\infty} \frac{x^2}{e^{k}+e^{x^2}}\,dx=-\frac{\sqrt{\pi }}{2} e^{-k} \,\text{Li}_{\frac{3}{2}}\left(-e^k\right)$$

I would like to add a different approach, directly connecting the integral to the field of polylogarithms (and I don't know why it wasn't mentioned within the two excellent answers already given).

Clearly, the integrand is an even function of $x$ and thus we only need to be concerned about the interval $[0,\infty]$. Enforcing the substitution $x^2\mapsto x$ reveals the integral representation of the Polylogarithm Function as \begin{align*} \int_{-\infty}^\infty \frac{x^2}{e^{x^2}+e^k}\mathrm dx&=2\int_0^\infty \frac{x^2}{e^{x^2}+e^k}\mathrm dx\\ &=\int_0^\infty\frac{x^\frac12}{e^x+e^k}\mathrm dx\\ &=-e^{-k}\int_0^\infty\frac{x^{\frac32-1}}{e^x/(-e^k)-1}\mathrm dx\\ &=-e^{-k}\left[\Gamma\left(\frac32\right)\operatorname{Li}_\frac32\left(-e^k\right)\right] \end{align*}

$$\therefore~\int_{-\infty}^\infty \frac{x^2}{e^{x^2}+e^k}\mathrm dx~=~-\frac{\sqrt\pi }2e^{-k}\operatorname{Li}_\frac32\left(-e^k\right)$$

The last step follows by basic properties of the Gamma Function. Anyway, note that the conditions - under which this representation holds - are met and at the same time the final expression matches the ones given by the other answers. Integral representation, as the one used here, are quite helpful when dealing with integrals of the "Mellin type" as they are sometimes called (due to their relation to the Mellin Transform).