I am not aware of any general results, but in this paper the authors have some thoughts for Gaussian variational approximations (GVAs) for generalized linear mixed model (GLMMs). Let $\vec y$ be the observed outcomes, $X$ be a fixed effect design matrix, $Z$ be a random effect design, denote an unknown random effect $\vec U$, and consider a GLMM with densities:
$$
\begin{align*}
f_{\vec Y\mid\vec U} (\vec y;\vec u) &=
\exp\left(\vec y^\top(X\vec\beta + Z\vec u)
- \vec 1^\top b(X\vec\beta + Z\vec u)
+ \vec 1^\top c(\vec y)\right) \\
f_{\vec U}(\vec u) &= \phi^{(K)}(\vec u;\vec 0, \Sigma) \\
f(\vec y,\vec u) &= f_{\vec Y\mid\vec U} (\vec y;\vec u)f_{\vec U}(\vec u)
\end{align*}
$$
where I use the same notation as in the paper and $\phi^{(K)}$ is a $K$-dimensional multivariate normal distribution density function.
Using a Laplace Approximation
Let
$$
g(\vec u) = \log f(\vec y,\vec u).
$$
Then we use the approximation
$$
\log\int \exp(g(\vec u)) d\vec u \approx
\frac K2\log{2\pi - \frac 12\log\lvert-g''(\widehat u)\rvert}
+ g(\widehat u)
$$
where
$$
\widehat u = \text{argmax}_{\vec u} g(\vec u).
$$
Using a Gaussian Variational Approximation
The lower bound in the GVA with a mean
$\vec\mu$ and covariance matrix $\Lambda$ is:
$$
\begin{align*}
\int \exp(g(\vec u)) d\vec u &\approx
\vec y^\top(X\vec\beta + Z\vec\mu)
- \vec 1^\top B(X\vec\beta + Z\vec\mu, \text{diag}(Z\Lambda Z^\top)) \\
&\hspace{25pt}+ \vec 1^\top c(\vec y) + \frac 12 \Big(
\log\lvert\Sigma^{-1}\rvert + \log\lvert\Lambda\rvert
-\vec\mu^\top\Sigma^{-1}\vec\mu \\
&\hspace{25pt} - \text{trace}(\Sigma^{-1}\Lambda)
+ K \Big) \\
B(\mu,\sigma^2) &= \int b(\sigma x + \mu)\phi(x) d x
\end{align*}
$$
where $\text{diag}(\cdot)$ returns a diagonal matrix.
Comparing the Two
Suppose that we can show that $\Lambda\rightarrow 0$ (the estimated conditional covariance matrix of the random effects tends towards zero). Then the lower bound (disregarding a determinant) tends towards:
$$
\begin{align*}
\int \exp(g(\vec u)) d\vec u &\approx
\vec y^\top(X\vec\beta + Z\vec\mu)
- \vec 1^\top b(X\vec\beta + Z\vec\mu) \\
&\hspace{25pt}+ \vec 1^\top c(\vec y) + \frac 12 \Big(
\log\lvert\Sigma^{-1}\rvert
-\vec\mu^\top\Sigma^{-1}\vec\mu + K\Big) \\
&= g(\vec\mu) + \dots
\end{align*}
$$
where the dots do not depend on the model parameters, $\vec\beta$ and $\Sigma$. Thus, maximizing over $\vec\mu$ yields $\vec\mu\rightarrow \widehat u$. Then the only difference between the Laplace approximation and the GVA is a
$$
- \frac 12\log\lvert -g''(\widehat u)\rvert
$$
term. We have that
$$
-g''(\widehat u) = \Sigma^{-1} + Z^\top b''(X\vec\beta + Z\vec u)Z
$$
where the derivatives are with respect to $\vec\eta = X\vec\beta + Z\vec u$. This does not tend towards zero as the conditional distribution of the random effects becomes more peaked. However, still very hand wavy, it may cancel out with the
$$
\frac 12\log\lvert\Lambda\rvert = -\frac 12\log\lvert\Lambda^{-1}\rvert
$$
term we disregarded in the lower bound. The first order condition for $\Lambda$ is:
$$
\Lambda^{-1} = \Sigma^{-1} + Z^\top B^{(2)}(X\vec\beta + Z\vec\mu, \text{diag}(Z\Lambda Z^\top)Z
$$
where
$$
B^{(2)}(\mu,\sigma^2) = \int b''(\sigma x+ \mu)\phi(x) dx.
$$
Thus, if $\vec\mu \approx \widehat u$ and $\Lambda \approx 0$ then:
$$
\Lambda^{-1} \approx \Sigma^{-1} + Z^\top b''(X\vec\beta + Z\vec u)Z
$$
and the Laplace approximation and the GVA yield the same approximation of the log marginal likelihood.
Notes
Do also see the annals paper Ryan Warnick mentions.