Your probability measure is a product measure, so by $$\text{Var}_y(\langle z,X\rangle) = \sum_{i=1}^nz_i^2\text{Var}_{y_i}(X_i)$$ everything reduces to the 1d case. Let $q_y(dx)=e^{xy-\varphi(x)-C(y)}dx$ be one of the marginals, where $y\in\mathbb{R}$ and $C(y)$ is chosen such that $q_y$ is normalized, and denote by $U_y$ the 1d r.v. with distributon $q_y$. It can be shown that $C(y)-C(0)$ is the cumulant-generating function of $U_0$, and $\text{Var}_y(U_y)=C''(y)$. Thus the variance has a local maximum at $y=0$ iff the third cumulant $\kappa_0^{(3)}$ of $U_0$ vanishes and the fourth cumulant $\kappa_0^{(4)}$ is negative. Indeed, if we denote by $m_y^{(k)}$ the $k$'th moment of $U_y$, we have $$\kappa_0^{(3)}=m_0^{(3)}-3m_0^{(1)}m_0^{(2)}+2\left(m_0^{(1)}\right)^3=0,$$ because $p_0$ is symmetric and therefore the first and third moments are zero. Similarly, dropping the odd-numbered moments, we have $$\kappa_0^{(4)}=m_0^{(4)}-3\left(m_0^{(2)}\right)^2.$$ Therefore assuming that $\varphi(x)>\lambda x^2$ for some $\lambda>0$ $(\star)$, then $$\kappa_0^{(4)}=\int_{\mathbb{R}^2}\left(x^4-3x^2y^2\right)e^{-\varphi(x)-\varphi(y)-2C(0)}dx dy < e^{-2C(0)}\int_{\mathbb{R}^2}\left(x^4-3x^2y^2\right)e^{-\lambda x^2-\lambda y^2}dx dy=0.$$ In order to show that this is the global maximum, let's show that $C''(y)$ is concave, i.e. $C^{(4)}(y)<0$ for all $y\in\mathbb{R}$. In fact, since $$C(y)=\log\int_{-\infty}^{\infty}e^{xy-\varphi(x)}dx,$$ $C^{(4)}(y)$ is the 4'th cumulant of $U_y$, i.e. $$C^{(4)}(y)\propto\int_{\mathbb{R}^4}\left(x_1^4-4x_1^3x_2-3x_1^2x_2^2+12x_1^2x_2x_3-6x_1x_2x_3x_4\right)e^{\sum_{i=1}^4(x_iy-\varphi(x_i))}dx\\<\int_{\mathbb{R}^4}\left(x_1^4-4x_1^3x_2-3x_1^2x_2^2+12x_1^2x_2x_3-6x_1x_2x_3x_4\right)e^{\sum_{i=1}^4(x_iy-\lambda x_i^2)}dx=0.$$
$(\star)$ This is slightly different than your 'increasing convexity' assumption, but seems sufficiently close.