I think the answer goes along these lines: Essentially it is Fubini and Substitution for multivariate integrals.
Let us assume that $d=2, \Sigma=\text{Id}$ and $\mu = 0$. Notice that the density takes the shape
$$\mathcal{N}(\Theta, \text{Id}, 0) = \text{const.}~ e^{-|\Theta|^2/2}$$
in this case. I will ignore the constant to keep the formulae simple. Take a rotation matrix
$$ A = \begin{pmatrix} \cos(\alpha) & -\sin(\alpha) \\ \sin(\alpha) & \cos(\alpha) \end{pmatrix}$$
such that $y_i A^T x_i = (c,0)^T$ where $c \neq 0$. If we execute the substitution (see https://en.wikipedia.org/wiki/Integration_by_substitution#Substitution_for_multiple_variables) $\Theta \mapsto A\Theta$ then we get that
$$\int_{\mathbb{R}^2} \phi(y_ix_i^T\Theta) \mathcal{N}(\Theta, \text{Id}, 0) d\Theta = \int_{\mathbb{R}^2} \phi(y_ix_i^TA\Theta) e^{-|A\Theta|^2/2} d\Theta$$
but notice that a vector $\Theta$ and its rotated variant $A\Theta$ have precisely the same length, i.e. $|A\Theta| = |\Theta|$ so that this simplifies to
$$\int_{\mathbb{R}^2} \phi((y_iAx_i)^T\Theta) e^{-|\Theta|^2/2} d\Theta = \int_{\mathbb{R}^2} \phi(c\Theta_1) e^{-\Theta_1^2/2} e^{-\Theta_2^2/2} d\Theta$$
At this point we use the Theorem of Fubini (https://en.wikipedia.org/wiki/Fubini%27s_theorem) which allows us to write the integral over the total domain $\mathbb{R}^2$ of $\Theta$ as an iterated integral
$$\int_{\mathbb{R}} \int_{\mathbb{R}} \phi(c\Theta_1) e^{-\Theta_1^2/2} e^{-\Theta_2^2/2} d\Theta_2 d\Theta_1 = \left( \int_{\mathbb{R}} e^{-\Theta_2^2/2} d\Theta_2 \right ) \cdot \left( \int_{\mathbb{R}} \phi(c\Theta_1) e^{-\Theta_1^2/2} d\Theta_1\right)$$
Now if we do not forget about the normalizing constant above and drag it along all the time then we have $(2\pi)^{-1}$ in front which we write as $\sqrt{2\pi}^{-1} \sqrt{2\pi}^{-1}$. The expression becomes
$$\left( \sqrt{2\pi}^{-1} \int_{\mathbb{R}} e^{-\Theta_2^2/2} d\Theta_2 \right ) \cdot \left( \sqrt{2\pi}^{-1} \int_{\mathbb{R}} \phi(c\Theta_1) e^{-\Theta_1^2/2} d\Theta_1\right)$$
The first integral vanishes (because it integrates to the constant $\sqrt{2\pi}$) and the second integral is what you wanted: an integral over the logistic function "in a univariate Gaussian expectation". If the amount of dimensions is more than $2$ then everything still works but gets a little messier. If $\Sigma$ is not the identity matrix and/or $\mu \neq 0$ then $\Sigma$ can be diagonalized unitarily (basic linear algebra) so you use some substitution $\Theta \mapsto U\Theta - \mu$ or so in advance in order to force the situation $\Sigma=\text{Id}$ and $\mu=0$ just as you do in the case of a univariate normal.