Multivariate Gaussian distribution is usually employed to decode real-valued data and Bernoulli distribution is usually employed to decode binary data in VAEs as claimed in your own reference.

${\displaystyle \ln p_{\theta }(x|z)}$ is implemented as ${\displaystyle -{\frac {1}{2}}\|x-D_{\theta }(z)\|_{2}^{2}}$, since that is, up to an additive constant, what ${\displaystyle x\sim {\mathcal {N}}(D_{\theta }(z),I)}$ yields. That is, we model the distribution of ${\displaystyle x}$ conditional on ${\displaystyle z}$ to be a Gaussian distribution centered on ${\displaystyle D_{\theta }(z)}$. The distribution of ${\displaystyle q_{\phi }(z|x)}$ and ${\displaystyle p_{\theta }(z)}$ are often also chosen to be Gaussians

In summary Gaussian, though may be limiting in modeling complex data such as GMMs, is often used to derive closed-form loss functions or solutions for real-valued data cases.