So, after a little while, I am able to answer my own question :) According to the definition of conditional density, we have
$$ f(\vec{x_a} | \vec{x_b}) = \frac{f(\vec{x_a}, \vec{x_b})}{f(\vec{x_b})} = \sum_{k=1}^{n} \frac{\pi_k}{f(\vec{x_b})} \mathcal{N}(\vec{x} | \mu_k, \Sigma_k). $$
where $ f(\vec{x_b}) $ is the marginal distribution. This can be calculated in terms of 1) the marginals and conditionals of multivariate Gaussians 2) the marginals of multivariate Gaussian mixtures.
Suppose that $ \mu $ and $ \Sigma $ is partitioned as
$$ \mu = (\mu_a, \mu_b) \text{ and } \Sigma = \begin{pmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ba} & \Sigma_{bb} \end{pmatrix}. $$
According to Bishop, Pattern Recognition and Machine Learning, Chapter 2.3.1-2.3.2, the conditional distributions of multivariate Gaussians can be calculated as
$$ \mathcal{N}(\vec{x_a} | \vec{x_b}, \mu, \Sigma) = \mathcal{N}(\vec{x_a} | \mu_{a|b}, \Sigma_{a|b}), $$
where
$$
\begin{aligned}
\mu_{a|b} & = \mu_a + \Sigma_{ab}\Sigma_{bb}^{-1}(x_b - \mu_b) \\
\Sigma_{a|b} & = \Sigma_{aa} - \Sigma_{ab}\Sigma_{bb}^{-1} \Sigma_{ba},
\end{aligned} $$
while the marginal is simply given by $ \mathcal{N}(\vec{x_b} | \mu_b, \Sigma_{bb}) $. Regarding the marginal of the MGM, we simply have
$$ f(\vec{x_b}) = \sum_{k=1}^{n} \pi_k \mathcal{N}(\vec{x_b} | \mu_{k, b}, \Sigma_{k, bb}), $$
so by continuing the calculation, we have
$$
\begin{aligned}
f(\vec{x_a} | \vec{x_b}) &= \sum_{k=1}^{n} \frac{\pi_k}{f(\vec{x_b})} \mathcal{N}(\vec{x} | \mu_k, \Sigma_k) \\
& = \sum_{k=1}^{n} \frac{\pi_k \mathcal{N}(\vec{x_b} | \mu_{k, b}, \Sigma_{k, bb})}{f(\vec{x_b})} \mathcal{N}(\vec{x_a} | \mu_{k, a|b}, \Sigma_{k, a|b}) \\
& = \sum_{k=1}^{n} \frac{\pi_k \mathcal{N}(\vec{x_b} | \mu_{k, b}, \Sigma_{k, bb})}{\sum_{l=1}^{n} \pi_l \mathcal{N}(\vec{x_b} | \mu_{l, b}, \Sigma_{l, bb})} \mathcal{N}(\vec{x_a} | \mu_{k, a|b}, \Sigma_{k, a|b}).
\end{aligned}
$$
Thus, it is indeed a multivariate Gaussian mixture with components $ \mathcal{N}(\vec{x_a} | \mu_{k, a|b}, \Sigma_{k, a|b}) $ and mixing coefficients $$ \pi_{k, a|b} = \frac{\pi_k \mathcal{N}(\vec{x_b} | \mu_{k, b}, \Sigma_{k, bb})}{\sum_{l=1}^{n} \pi_l \mathcal{N}(\vec{x_b} | \mu_{l, b}, \Sigma_{l, bb})}. $$