Assume that $X_1$ and $X_2$ are random variables with Gaussian joint distribution
$$ \begin{bmatrix} X_1 \\ X_2 \end{bmatrix} \sim \mathcal{N}( \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}, \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix} ). $$
Furthermore, Assume that $X_2$ and $X_3$ have Gaussian joint distribution
$$ \begin{bmatrix} X_2 \\ X_3 \end{bmatrix} \sim \mathcal{N}( \begin{bmatrix} \mu_2 \\ \mu_3 \end{bmatrix}, \begin{bmatrix} \Sigma_{22} & \Sigma_{23} \\ \Sigma_{32} & \Sigma_{33} \end{bmatrix} ). $$
To my understanding, this does not imply that $[X_1 \ X_2 \ X_3]^T$ has a multivariate Gaussian joint distribution and also does not imply anything about ${\rm cov}(X_1,X_3)$.
Now consider the following method of sampling:
First, let $x_3$ be sampled from $X_3 \sim \mathcal{N}(\mu_3,\Sigma_{33})$. Then, let $x_2$ be sampled from the conditional distribution $(X_2 | X_3 = x_3) \sim \mathcal{N}(\mu_{2 | 3},\Sigma_{2 | 3})$, where
$$ \mu_{2 | 3} = \mu_2 + \Sigma_{23} \Sigma_{33}^{-1} (x_3 - \mu_3), \\ \Sigma_{2 | 3} = \Sigma_{22} - \Sigma_{23} \Sigma_{33}^{-1} \Sigma_{32}, $$
which are obtained from the conditional distribution of multivariate Gaussians. If I'm not mistaken, this is equivalent to $[x_1, x_2]^T$ being sampled from $\mathcal{N}( \begin{bmatrix} \mu_2 \\ \mu_3 \end{bmatrix}, \begin{bmatrix} \Sigma_{22} & \Sigma_{23} \\ \Sigma_{32} & \Sigma_{33} \end{bmatrix} )$.
Finally, I sample $x_1$ from the conditional distribution $(X_1 | X_2 = x_2) \sim \mathcal{N}(\mu_{1 | 2},\Sigma_{1 | 2})$, where
$$ \mu_{1 | 2} = \mu_1 + \Sigma_{12} \Sigma_{22}^{-1} (x_2 - \mu_2), \\ \Sigma_{2 | 3} = \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1} \Sigma_{32}. $$
My question is what distribution is $[x_1,x_2,x_3]^T$ being sampled from? Is it from a multivariate Gaussian of the form $\mathcal{N}( \begin{bmatrix} \mu_1 \\ \mu_2 \\ \mu_3 \end{bmatrix}, \begin{bmatrix} \Sigma_{11} & \Sigma_{12} & {\rm cov}(X_1,X_3) \\ \Sigma_{21} & \Sigma_{22} & \Sigma_{23} \\ {\rm cov}(X_3,X_1) & \Sigma_{32} & \Sigma_{33} \end{bmatrix})$? If so, is there some inherent assumption in the method of sampling that defines what ${\rm cov}(X_1,X_3)$ is?
For example, here's $10^5$ samples with $\mu_1 = \mu_2 = \mu_3 = 0$, $\Sigma_{11} = \Sigma_{22} = \Sigma_{33} = 1$, and $\Sigma_{12} = \Sigma_{21} = \Sigma_{23} = \Sigma_{32} = 0.9$:
