Global minimum at $0$ for $p<2$ is trivial in the sense that you can see it without writing a single equation or inequality (i.e., a non-trivial formula with two sides and a sign in between comparing them).
Observation 1 $e^{-|x|^p}$ is some weighted average of $e^{-ax^2}$ with positive $a$.
Observation 2 The thingy you are interested in is just the variance of $x$ with respect to the probability measure $\mu_y$ whose density is proportional to $e^{-x^2}e^{-|x+y|^p}$. By observation 1, this measure is a mixture of the probability measures $\mu_{a,y}$ with densities proportional to $e^{-x^2}e^{-a(x+y)^2}$. The weight of $\mu_{a,y}$ in that mixture is proportional to something independent of $y$ times $e^{-\frac a{a+1}y^2}$, i.e., when we move $y$ away from the origin, the measures $\mu_{a,y}$ with lower $a$ gain more weight in the composition.
Observation 3. The variance of $x$ with respect to $\mu_{a,y}$ is independent of $y$ and decreases in $a$. To be exact, it is just inversely proportional to $1+a$.
Observation 4. The variance in the mixture is at least the mixture of the variances, which is minimized at $0$ by the independence of individual variances of $y$, observation 3, and the last sentence of observation 2. Also at $y=0$ we have equality because all means are at $0$ by symmetry.
The end.
I wish I could come up with an equally simple argument for $p>2$, but, alas, I don't have one at the moment.
Edit OK, I guess I finally figured it out. We'll prove the following statement. Let $\varphi$ be an even convex function such that $\varphi''$ increases on $[0,+\infty)$. Let $p_y(x)$ be the probability density proportional to $e^{yx}e^{-\varphi(x)}$. Then the variance of $x$ with respect to the corresponding probability measure is a non-increasing function of $y$ for $y>0$.
Indeed, let's differentiate in $y$. We have $p_{y+\delta y}(x)$ proportional to $p_y(x)(1+x\delta y)$ (the linearization of $e^{x\delta y}$), but if we leave it at that, the mass will change to $1+\delta y\int p_y(x)xdx=1+c\delta y$ where $c$ is the expectation of $x$ with respect to $p_y$, so we need to compensate by dividing by that factor, which will result in the linearization
$$
p_{y+\delta y}(x)=(1+(x-c)\delta y)p_y(x)
$$
Taking the linearization of the variance of $x-c$ now (which is the same as the variance of $x$ but easier to compute), we see that what we need to show is that
$$
\int (x-c)^3p_y(x)dx\le 0.
$$
We will show even that $\int_{x:|x-c|>a}(x-c)p_y(x)\le 0$ for all $a\ge 0$. The equality obviously holds for $a=0$ (the definition of $c$) and for $a=+\infty$. The derivative in $a$ is just $a(p_y(c-a)-p_y(c+a))$. I claim now that it can change sign only once for $a>0$ and that change is from $-$ to $+$. That is equivalent to saying that $\Phi(a)=\varphi_y(c-a)-\varphi_y(c+a)$, where $\varphi_y(x)=\varphi(x)-yx$, can change sign only once for $a>0$ and that change is from $+$ to $-$.
We clearly have $c>0$ for $y>0$. Therefore, the point $c-a$ is always closer to the origin than $c+a$ for $a>0$ whence, by the assumed property of the second derivative of $\varphi$ (the linear term coming from $yx$ can do nothing with the second derivative), $\Phi''<0$. Thus, if we start from $\Phi(0)=0$ in the positive direction, we change sign once from $+$ to $-$ as promised. If we had started in the negative direction, we would never be able to change sign at all and the integral we are interested in would be monotone all the way, which is ridiculous, because its values at $0$ and $+\infty$ are both $0$, so that case is impossible.
That takes care of $p>2$. If the second derivative of $\varphi$ were decreasing, then all inequalities would be reversed, so we can cover $p<2$ by this method as well and even get unimodality you see on the pictures, not just the global minimum at $0$.
I hope I haven't made any stupid mistake in the computations but, since it is nearly midnight here now, you'd better check them carefully :-)
